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ADVANCED  COMPUTATIONAL  TECHNIQUES 


IN  REGIONAL  WAVE  STUDIES 


DISCUSSION 

The  purpose  of  this  contract  was  to  perform  computational  studies  that  would 
increase  our  understanding  of  the  Lg  arrival  and  its  coda,  especially  the  use  of  the 
coda  to  define  a  path  specific  attenuation  operator  for  the  Lg  arrival.  This  is  a  very 
difficult  problem  which  requires  the  incorporation  of  realistic  scattering  mechanisms 
for  forward  modeling  of  the  phenomena. 

Work  done 

It  was  realized  that  stable,  well  understood,  synthetic  seismograms  for  simpler 
homogeneous  plane  layers  is  required  first.  Toward  this  end,  programs  developed  for 
other  projects  were  documented  and  distributed  through  the  Center  for  Seismic  Studies 
for  use  by  other  researchers.  These  programs  consisted  of  over  80,000  lines  of  code 
and  800  pages  of  documentation.  Much  effort  was  spent  refining  this  package,  espe¬ 
cially  those  portions  written  by  students,  to  ensure  legibility  and  numerically  stable 
results. 

The  set  of  six  volumes  of  "Computer  Programs  in  Seismology"  were  distributed 
by  the  Center  for  Seismic  Studies  to  the  following  organizations: 


Organization 

Contact 

ENSCO  (FLA) 

H.  Ghalib 

ENSCO  (VA) 

Z.  Der 

SMU 

B.  Stump 

U.  Ruhr 

H.-P.  Haijes 

USGS 

S.  Sipkin 

S.  Radiomana 

B.  Massinon 

AFTAC 

N.  Yacoub 

In  addition,  other  DARPA/AFGL  contractors  have  obtained  copies  directly  from  Saint 
Louis  University. 

Toward  the  latter  part  of  the  project,  the  computer  programs  were  stable  enough 
to  attempt  to  look  at  real  data.  As  reported  in  Quarterly  Management  Report  No.  8, 
the  Lg  and  surface-wave  signals  from  a  strip  mining  blast  250  km  away  in  western 
Kentucky  were  used  to  infer  a  source  yield  of  about  0.1  kT.  The  implication  of  this 
result  is  that  more  than  one  wavetype  at  regional  distances  can  be  used  to  obtain  a 
robust  seismic  estimate  yield  of  small  explosions. 

One  of  the  unknowns  is  the  value  of  Q  p  in  the  upper  3  km  of  the  earth.  We  have 
been  looking  at  data  from  refraction  profiles  in  Maine,  from  shallow  seismic  investiga¬ 
tions  in  flood  plains,  and  from  regionally  recorded  mining  blasts  to  define  both  the 
robust  analysis  procedures  required  and  the  shallow  shear-wave  velocity  and  Q  p. 


Three  dissertations  were  completed: 

Russell,  D.  R.,  (1987).  Multi-channel  processing  of  dispersed  surface  waves. 

Shieh,  C.-F.,  (1988).  Polarization  analysis  of  complex  seismic  wave  field. 

M.  L.  lost  (1989).  Long  period  strong  ground  motions  and  response  spectra  of  large 
historical  earthquakes  in  the  central  and  eastern  United  States  from  kinematic  source 
models,  Ph.  D.  Dissertation,  Saint  Louis  University. 

David  Russell  is  currently  with  AFTAC,  performing  work  related  to  that  sup¬ 
ported  by  DARPA/AFGL.  Michael  Jost  has  taken  a  position  at  the  University  of  the 
Ruhr,  working  with  the  new  GERES  S  data. 

The  dissertation  work  emphasized  the  development  and  use  of  advanced  computa¬ 
tional  techniques  for  studying  regional  seismic  phases.  Russell  (1987)  implemented 
surface-wave  inversion,  single  and  multi-channel  phase  matched  filtering  for  the 
analysis  of  dispersed  surface  waves.  Shieh  (1988)  focused  on  polarization  filtering, 
with  emphasis  on  exploration  data,  but  he  also  considered  representative  teleseismic 
signals.  Jost  (1989)  made  use  of  the  programs  to  estimate  low  frequency  ground 
motion  at  short  distances  from  large  earthquakes.  He  found  it  necessary  to  modify  a 
scaling  relationship  for  continental/plate  interiors  proposed  by  Nuttli  (Jost’s  paper  is 
attached).  His  work  provides  a  better  definition  of  how  intra-continental  earthquakes 
should  behave. 

Recommendations 

The  discrimination  and  quantification  problems  become  difficult  for  small  events. 
On  the  other  hand,  the  possibility  of  new  data  sources  at  regional  distances  permits 
using  previously  ignored  signals.  Unfortunately,  these  regional  signals  will  travel  large 
distances  through  a  heterogeneous  crust.  Attempts  must  be  made  to  use  the  data  to 
define  path  dependent  attenuation  operators  and  to  understand  the  effects  of  both 
source  and  path  heterogeneity  on  the  observed  signal.  Even  though  moment  tensor 
inversion  may  not  be  used  too  much  for  large  explosions,  such  inversion  of  regional 
phase  data,  especially  surface-wave,  may  provide  added  constraints  on  whether  the 
small  explosive  source  is  point  or  distributed. 
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A  Student's  Guide  to  and  Review  of  Moment  Tensors 

M.  L.  Jost  and  R.  B.  Herrmann 

Department  of  Earth  and  Atmospheric  Sciences 
Saint  Louis  University 
P  O  Box  8009 
St.  Louis,  A/O  6S156 

ABSTRACT 

A  review  of  a  moment  tensor  for  describing  a  general  seismic  point  source  is 
presented  to  show  a  second  order  moment  tensor  can  be  related  to  simpler  seismic 
source  descriptions  such  as  centers  of  expansion  and  dojble  couples.  A  review  of 
literature  is  followed  by  detailed  algebraic  expansions  of  the  moment  tensor  into 
isotropic  and  deviatoric  components.  Specific  numerical  examples  are  provided  in 
the  appendices  for  use  in  testing  algorithms  for  moment  tensor  decomposition. 


INTRODUCTION 

A  major  research  interest  in  seismology  is  the 
description  of  the  physics  of  seismic  sources.  A  common 
approach  is  the  approximation  of  seismic  sources  by  a 
model  of  equivalent  forces  that  correspond  to  the  linear 
wave  equations  neglecting  non-linear  effects  in  the  near 
source  region  (Geller,  1076;  Aki  and  Richards,  1080;  Ken- 
nett,  1083;  Bullen  and  Bolt,  1085).  Equivalent  forces  are 
defined  as  producing  displacements  at  the  earth's  surface 
that  are  identical  to  those  from  the  actual  forces  of  the 
physical  process  at  the  source.  The  equivalent  forces  are 
determined  from  observed  seismograms  that  contain  infor¬ 
mations  about  the  source  and  path  and  distortions  due  to 
the  recording.  Hence,  the  principle  problem  of  source  stu¬ 
dies  is  the  isolation  of  the  source  effect  by  correcting  for 
instrument  and  path. 

The  classical  method  of  describing  seismic  sources, 
having  small  dimensions  compared  to  the  wavelengths  of 
interest  (point  source  approximation)  is  by  their  strength 
(magnitudes,  seismic  moment)  and  their  fault  plane  solu¬ 
tion  (Honda,  1062;  Hirasawa  and  Stauder,  1065; 
Herrmann,  1975).  Recently,  seismic  moment  tensors  have 
been  used  routinely  for  describing  seismic  point  sources 
(e.g.  Kanamori  and  Given,  1082;  Dziewonski  and  Wood- 
house,  1083b;  Dziewonski  et  of.  ,  1083a-c,  108-la-c;  Giar- 
dini,  1084;  Ekstrom  and  Dziewonski,  1085;  Dziewonski 
et  at  ,  1085a-d,  1086a-c,  1087a-f;  Ekstrom  et  at. ,  1087; 
Sipkin,  1087;  PDE  monthly  listings  published  by  NEIS). 
Gilbert  (1070)  introduced  moment  tensors  for  calculating 
the  displacement  at  the  free  surface  which  can  be 
expressed  as  a  sum  of  moment  tensor  elements  times  the 
corresponding  Green’s  function.  An  elastodynamic 
Green’s  function  is  a  displacement  field  due  to  an  uni¬ 
directional  unit  impulse,  i.e.  the  Green's  function  is  the 
impulse  response  of  the  medium  between  source  and 
receiver.  The  response  of  the  medium  to  any  other  time 
function  is  the  convolution  (Arfken,  1085)  of  that  time 
function  with  the  impulse  response.  The  Green’s  function 
depends  on  source  and  receiver  coordinates,  the  earth 
model,  and  is  a  tensor  (Aki  and  Richards,  1080).  The 
linearity  between  the  moment  tensor  and  Green’s  function 
elements  was  first  used  by  Gilbert  (1073)  for  calculating 
moment  tensor  elements  from  observations  (moment  ten¬ 
sor  inversion).  The  concept  of  seismic  moment  tensors 


was  further  extended  by  Backus  and  Klulcahy  (1076),  and 
Backus  (1077a,  b).  Moment  tensors  can  be  determined 
from  free  oscillations  of  the  earth  (e.g.  Gilbert  and 
Dziewonski,  1975),  long-period  surface  waves  (e.g. 
McCowan,  1076;  Mendiguren,  1077;  Patton  and  Aki, 
1070;  Patton,  1080;  Kanamori  and  Given,  1981,  1082; 
Romanowicz,  1981;  Lay  et  at.,  1982;  Nakanishi  and 
Kanamori,  1082,  1084)  or  long-period  body  waves  (e.g. 
Stump  and  Johnson,  1077;  Strelitz,  1078,  1980;  Ward, 
1080a,  b;  Fitch  et  at,  1080;  Fitch,  1981;  Langston,  1081; 
Dziewonski  et  at,  1081;  Dziewonski  and  Woodhouse, 
1083a,  b).  Throughout  this  Student's  Guide,  we  will 
focus  on  second-rank,  time  independent  moment  tensors 
(Appendix  I).  We  refer  to  Dziewonski  and  Gilbert  (1974), 
Gilbert  and  Dziewonski  (1975),  Backus  and  Mulcahv 
(1076),  Backus  (1077a),  Stump  and  Johnson  (1977),  Strel¬ 
itz  (1080),  Sipkin  (1082),  and  Vasco  and  Jo'  nson  (1988) 
for  a  description  of  time  dependent  moment  tensors. 
Higher  order  moment  tensors  are  discussed  by  Backus  and 
Mulcahy  (1076),  Backus  (1077a,  b),  and  Dziewonski  and 
Woodhouse  (1083a). 

The  reason  that  moment  tensors  are  important  is 
that  they  completely  describe  in  a  first  order  approxima¬ 
tion  the  equivalent  forces  of  general  seismic  point  sources. 
The  equivalent  forces  can  be  correlated  to  physical  source 
models  such  as  sudden  relative  displacement  at  a  fault 
surface  (elastic  rebound  model  by  H.  F.  Reid,  1010), 
rapidly  propagating  metastable  phase  transitions  (Evison, 
1963),  sudden  volume  collapse  due  to  phase  transitions,  or 
sudden  volume  increase  due  to  explosions  (Kennctt,  1083; 
Va sco  and  Johnson,  1988)  The  equivalent  forces 
representing  a  sudden  displacement  on  a  fault  plane  form 
the  familiar  double  couple.  The  equivalent  forces  of  a  sud¬ 
den  change  in  shear  modulus  in  presence  of  axial  strain 
are  represented  by  a  linear  vector  dipole  (Knopoff  and 
Randall.  1970).  In  conclusion,  a  seismic  moment  tensor  is 
a  general  concept,  describing  a  variety  of  seismic  source 
models,  the  shear  dislocation  (double  couple  source)  being 
just  one  of  them. 

The  equivalent  forces  can  be  determined  from  an 
analysis  of  the  eigenvalues  and  eigenvectors  of  the 
moment  tensor  (Appendix  I).  The  sum  of  the  eigenvalues 
of  the  moment  tensor  describee  the  volume  change  in  the 
source  (isotropic  component  of  the  moment  tensor).  If  the 
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sum  is  positive,  the  isotropic  component  is  due  to  an 
explosion.  The  source  has  an  implosive  component  if  the 
sum  is  negative.  If  the  sum  of  the  eigenvalues  vanishes, 
then  the  moment  tensor  has  only  deviatoric  components. 
The  deviatoric  moment  tensor  represents  a  pure  double 
couple  source  if  one  eigenvalue  equals  zero.  If  none  of  the 
eigenvalues  vanishes  and  their  sum  still  equals  zero,  the 
moment  tensor  can  be  decompcsed  into  a  major  and 
minor  double  couple  (Kanamori  and  Given,  1081),  or  a 
double  couple  and  a  compensated  linear  vector  dipole 
(CLVD)  (Knopoff  and  Randall,  1970).  A  CLVD  is  a 
dipole  that  is  corrected  for  the  effect  of  volume  change, 
describing  seismic  sources  which  have  no  volume  change, 
net  force,  or  net  moment.  In  general,  a  complete  moment 
tensor  can  be  the  superposition  of  an  isotropic  component 
and  three  vector  dipoles  (or  three  CLVD's,  or  three  double 
couples,  Ben-Menahem  and  Singh,  1981). 

This  Student’s  Guide  is  an  extension  of  "A  student's 
guide  to  the  use  of  P-  and  S-  wave  data  for  focal  mechan¬ 
ism  determination"  (Herrmann,  1975).  Hence,  emphasis  is 
given  illustrating  the  relations  between  classical  fault 
plane  solutions  and  seismic  moment  tensors.  Addressing 
general  seismic  point  sources,  we  provide  examples  of 
moment  tensor  decompositions  into  basic  equivalent 
source  representations,  as  contributions  of  dipoles  or  dou¬ 
ble  couple  sources.  Clarification  of  terms  such  as  major 
and  minor  double  couple  or  compensated  linear  vector 
dipole  is  provided.  Moment  tensor  inversion  schemes  are 
briefly  summarized.  In  the  appendices,  examples  of  the 
use  of  notation  by  different  authors  are  given  along  with 
some  numerical  results  which  are  useful  for  testing  com¬ 
puter  programs.  Furthermore,  the  formulation  of  the 
basic  Green's  functions  by  Herrmann  and  Wang  (1985)  is 
connected  to  a  simple  moment  tensor  inversion  scheme. 

GENERAL  ELASTODYNAMIC  SOURCE 

By  using  the  representation  theorem  for  seismic 
sources  (Aki  and  Richards,  1980),  the  observed  displace¬ 
ment  dn  at  an  arbitrary  position  x  at  the  time  t  due  to  a 
distribution  of  equivalent  body  force  densities,  fk ,  in  a 
source  region  is 

OO 

d„(x,t)~  /  /  Gnk(x,l\T,t)  fk{r,t)  dV(r)  dt  ,  (1) 

-oo  V 

where  Gnk  are  the  components  of  the  Green's  function 
containing  the  propagation  effects,  and  V  is  the  source 
volume  where  (k  are  non-zero.  We  assume  the  summation 
convention  for  repeated  indices  (Arfken,  1985).  The  sub¬ 
script  n  indicates  the  component  of  the  displacement. 
Throughout,  we  will  use  the  following  coordinate  system 
(Figure  1):  The  x-axis  points  towards  north,  the  y-axis 
towards  east,  and  the  z-axis  down  (this  system  is  right 
handed).  Then,  ex,  ey,  and  e,  are  the  unit  vectors 
towards  north,  east,  and  vertically  down,  respectively. 

By  assuming  that  the  Green's  functions  vary 
smoothly  within  the  source  volume  in  the  range  of 
moderate  frequencies,  the  Green's  functions  can  be 
expanded  into  a  Taylor  series  around  a  reference  point  to 
facilitate  the  spatial  integration  in  (1)  (Konnett,  1983; 
Arfken,  1985).  The  expansion  is  usually  done  around  the 
centroid  r  =  £.  The  physical  source  region  is  character¬ 
ized  by  the  existence  of  the  equivalent  forces.  These  forces 


North 


Fig.  1.  Definition  of  the  Cartesian  coordinates  (x,y.z). 
The  origin  is  at  the  epicenter.  Strike  is  measured 
clockwise  from  north,  dip  from  horizontal  down, 
and  slip  counterclockwise  from  horizontal,  u  and  i> 
are  the  slip  vector  and  fault  normal,  respectively 
(modified  after  Aki  and  Richards.  1980). 


arise  due  to  differences  between  the  model  stress  and  the 
actual  physical  stress  (stress  glut,  Backus  and  Mulcahy, 
1976).  Outside  the  source  region,  the  stress  glut  vanishes 
as  do  the  equivalent  forces.  The  centroid  of  the  stress 
glut  is  then  a  weighted  mean  position  of  the  physical 
source  region  (Backus,  1977a:  Aki  and  Richards.  1980; 
Dziewonski  and  Woodhouse,  1983a).  It  seems  that  the 
centroid  of  the  stress  glut  gives  a  better  position  for  the 
equivalent  point  source  of  an  earthquake  than  the  hypo- 
center  which  describes  just  the  position  of  rupture  initiali¬ 
zation.  The  Taylor  series  expansion  of  the  components  of 
the  Green's  function  around  this  new  reference  point  is 

Gnk(x,t;rJ)-  (2) 


The  comma  be*  ween  indices  in  (2)  describes  partial  deriva¬ 
tives  with  respect  to  the  coordinates  after  the  comma. 
We  define  the  components  of  the  time  dependent  force 
moment  tensor  as  ; 


A/*;,  >.(€.* )  “  /  (r, ,-*>,)  (r.O  dV  (3) 

v 

If  conservation  of  linear  momentum  applies,  such  as  for  a 
source  in  the  interior  of  a  body,  then  a  term  in  Mk  does 
not  exist  in  (3).  With  the  Taylor  expansion  (2)  and  the 
definition  of  the  time  dependent  moment  tensor  (3).  the 
displacement  (1)  can  be  written  as  a  sum  of  terms  which 
resolve  additional  details  of  the  source  (multipole  expan- 
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sion,  Backus  and  Mulcahy,  1976;  Stump  and  Johnson,  1082).  Neglecting  higher  order  terms,  we  get  (Stump  and 


1977;  Aki  and  Richards,  1080;  Kennett,  1983;  Dziewonski 
and  Woodhouse,  1083a;  Vasco  and  Johnson,  1088): 

(<) 

m- 1  m  • 

where  *  denotes  the  temporal  convolution.  By  using  a 
seismic  signal  that  has  much  longer  wavelengths  than  the 
dimensions  of  the  source  (point  source  approximation),  we 
need  to  consider  only  the  first  term  in  (4)  (Backus  and 
Mulcahy,  1976;  Stump  and  Johnson,  1977).  Note,  that 
single  forces  will  not  be  present  in  (4)  if  there  are  no 
externally  applied  forces  (indigenous  source).  The  tc’al 
force,  linear  and  angular  momentum  must  vanish  for  the 
equivalent  forces  of  an  indigenous  source  (Backus  and 
Mulcahy,  1976).  The  conservation  of  angular  momentum 
for  the  equivalent  forces  leads  to  the  symmetry  of  the 
seismic  moment  tensor  (Gilbert,  107C). 

We  assume  that  all  components  of  the  time  depen¬ 
dent  seismic  moment  tensor  in  (4)  have  the  same  time 
dependence  s(l)  (synchronous  source,  Silver  an<'  Jordan, 


z  z  Z 

Fig.  2.  The  nine  generalized  couples  representing  G„k  j  in  (51  (modified 
after  Aki  and  Richards,  1980). 
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Johnson, 1977) 

d„(x,t)  =  Mt, \  CHl  ]  *  *(0  )  (5) 

Mk)  are  constants  representing  the  components  of  the 
second  order  seismic  moment  tensor  M,  usually  termed 
Ihe  moment  tensor.  Note  that  the  displacement  dn  is  a 
linear  function  of  the  moment  tensor  elements  and  the 
terms  in  the  square  brackets.  If  the  source  time  function 
s(()  is  a  delta  function,  the  only  term  left  in  the  square 
brackets  is  Gnk  ,  describing  nine  generalized  couples.  The 
derivative  of  a  Green’s  function  component  with  respect 
to  the  souice  coordinate  is  equivalent  to  a  single  couple 
with  arm  in  the  direction.  For  k  =  j,  i.e.  force  in  the 
same  direction  as  the  arm,  the  generalized  couples  are  vec¬ 
tor  dipoles  (Figure  2;  Maruyama,  1964).  Thus,  the 
moment  tensor  component  Mk)  gives  the  excitation  of  the 
generalized  (k,j)  couple. 

DOUBLE  COUPLE  SOURCES 

The  moment  tensor  components  in  (5)  in  an  isotropic 
medium  for  a  double  couple  of  equivalent  forces  are  given 
by 
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Mkj  -  p  A  (  uk  Uj  +  u;-  i/k  )  ,  (6) 

where  p  is  the  shear  modulus,  A  is  the  area  of  the  fault 
plane,  u  denotes  the  slip  vector  on  the  fault  surface,  and 
v  is  the  vector  normal  to  the  fault  plane  (Aki  and 

Richards,  1980;  Ben-Menahem  and  Singh,  1081).  Note 
that  the  contributions  of  the  vector  of  the  fault  normal  u 
and  the  slip  vector  u  are  symmetric  in  (6).  From  the  sym¬ 
metry  of  M,  we  note  that  the  roles  of  the  vectors  u  and  v 
could  be  interchanged  without  affecting  the  displacement 
field;  i.e.  the  fault  normal  could  equivalently  be  the  slip 
vector  and  vice  versa.  This  well  known  far’  plane  -  auxi¬ 
liary  plane  ambiguity  cannot  be  resolved  from  the  seismic 
radiation  of  a  point  source.  Hence  studies  of  locations  of 
aftershocks,  surface  faulting,  rupture  directivity,  or  static 
final  displacements  (Backus,  1077a)  need  to  be  done  in 
order  to  resolve  this  ambiguity. 

The  term  vk  v.  4-  u}  uk  in  (6)  forms  a  tensor,  D, 
describing  a  double  couple.  This  tensor  is  real  and  sym¬ 
metric,  giving  real  eigenvalues  and  orthogonal  eigenvec¬ 
tors  (Appendix  1).  The  eigenvalues  are  proportional  to  (1, 
0,  -1).  Hence,  the  characteristic  properties  of  a  moment 
tensor  representing  a  double  couple  are  i)  one  eigenvalue 
of  the  moment  tensor  vanishes,  and  ii)  the  sum  of  the 
eigenvalues  vanishes,  i.e.  the  trace  of  the  moment  tensor  is 
zero  (t he  other  two  eigenvalues  are  constrained  to  equal 
magnitude  but  opposite  sign). 

Let  t,  b,  and  p  designate  the  orthogonal  eigenvectors 
to  the  above  eigenvalues  (Herrmann,  1975;  Backus,  1977a; 
Dziewonski  and  Woodhouse,  1983a). 


t  «=  (  v  -F  u  ) 

(7) 

b  =  v  X  u 

(8) 

p  =  ^  ^  ~  U  * 

(0) 

The  tensor  D  corresponding  to  the  terms  in  the  brackets 
in  (6)  can  be  diagonalized  (principal  axis  transformation, 
see  Appendix  I),  where  the  eigenvectors  give  the  directions 
of  the  principal  axes.  The  eigenvector  b  corresponding  to 
the  eigenvalue  zero  gives  the  null-axis,  the  eigenvector  t 
corresponding  to  the  positive  eigenvalue  gives  the  tension 
axis,  T,  and  the  eigenvector  p  corresponding  to  the  nega¬ 
tive  eigenvalue  gives  the  pressure  axis,  P,  of  the  tensor. 
These  axes  can  be  related  to  the  corresponding  axes  of  the 
fault  plane  solution,  since  we  are  focusing  on  pure  double 
couple  sources.  The  P-axis  is  in  the  direction  of  max¬ 
imum  compressive  motion  on  the  fault  surface;  the  T-axis 
is  the  direction  of  maximum  tensional  motion.  Note  that 
the  P-  and  T-axcs  inferred  from  the  motion  on  the  fault 
surface  are  not  necessarily  identical  to  the  axes  of 
maximum  tectonic  stress,  sine  the  motion  can  be  on  a 
preexisting  plane  of  weakness  rather  than  on  a  newly 
formed  fault  plane  that  would  correspond  to  the  max¬ 
imum  tectonic  stress  (McKenzie,  1969).  However,  this 
ambiguity  cannot  be  resolved  from  the  seismic  radiation. 
In  order  to  determine  the  direction  of  maximum  tectonic 
stress,  additional  geological  data  such  as  in  situ  stress 
measurements  and  frictional  forces  is  necessary.  1  ncking 
this  kind  of  information,  it  is  generally  assumed  that  the 
P-  and  T-  axes  found  from  the  seismic  wave  radiation  are 
somewhat  indicative  of  the  direction  of  tectonic  stress. 


The  double  couple  uk  t/j  +  u;  uk  can  equivalently  be 
described  by  its  eigenvectors  (Gilbert,  1973). 

uk  Vj  +  Uj  uk  -  tk  t}  -  pj  pk  (10) 

=  0.5  |(<t  +  pk)  (tj  -  Pj)  +  (<*  ~  Pk)  ((;  +  P,)\ 

Comparing  the  terms  in  (10),  we  find  the  relation  be'we-n 
tension  and  pressure  axes  and  slip  vector  and  fault  normal 


(Appendix  I): 

u  =  (  t  +  p  )  (11) 

‘,  =  W(‘~p)  (12) 

The  other  nodal  ri  me  is  defined  by 

u  —  (  t  -  p  )  (13) 

+P)  (H) 


If  strike,  4>,  dip,  5,  and  slip,  X,  of  the  faulting  are 
known,  the  slip  vector  u  and  the  fault  normal  u  are  given 
by  (Aki  and  Richards,  1980) 

u  =  u  (  cos  X  cos  4>  +  cos  5  sin  X  sin  )  e, 

■f  iT  (  cos  X  sin  —  cos  8  sin  X  cos  *)*,  (15) 

—  iT  sin  8  sin  X  e,  , 

where  iT  is  the  mean  displacement  on  the  fault  plane.  The 
fault  normal  v  is 

1/  =  —  sin 6  sin<l>  e,  +  sin5  cos<I>  ey  —  cos5  e,  (16) 

The  scalar  product  of  u  and  v  is  zero.  The  strike  of  the 
fault  plane,  4>,  is  measured  clockwise  from  north,  with  the 
fault  plane  dipping  to  the  right  when  looking  along  the 
strike  direction.  Equivalently,  the  hanging  wall  is  then  to 
the  right  (Figure  1).  The  dip,  8,  is  measured  down  from 
the  horizontal.  The  slip,  X,  is  the  angle  between  the  strike 
direction  and  the  direction  the  hanging  wall  moved  rela¬ 
tive  to  the  foot  wall  (the  slip  is  positive  when  measured 
counterclockwise  as  viewed  from  the  hanging  wall  side). 
The  range  of  the  fault  orientation  para  ieters  are 

0  <  <  2n,  0<5<y,  and  —  n  <  X  <  zr  (Herrmann, 

1975;  Aki  and  Richards,  1980).  The  scalar  seismic 
moment  is 

M„  *=  n  A  u  (17) 

Equation  (6)  together  with  (15),  (16),  and  (17)  lead  to 
the  Cartesian  components  of  the  symmetric  moment  ten¬ 
sor  in  terms  of  strike,  dip,  and  slip  angles. 

A/„  «■  — A/„( sin  8  cos  X  sin  2<fr  +  sin  2 8  sin  X  sin2  4>) 

Mfs  **  A/,, (sin  8  cos  X  sin  2$  —  sin  2 8  sin  X  cos2  <J>) 

A/„  =  A/0(sin  2 8  sin  X)  (18) 

«  A/„(sin  8  cos  X  cos  2$  +  0.5  sin  2 8  sin  X  sin  2$) 

A/„  —  —A/,  (cos  8  cos  X  cos  4>  +  c«.s  25  sin  X  sin  ♦) 

Afy>-  -  -M,  (cos  5  cos  X  s,n  <£  —  cos  25  sin  X  cos  $) 

Different  notation  of  the  moment  tensor  elements  are  dis¬ 
cussed  in  App.r.dix  II.  In  Appendix  III,  several  s:mp!e 
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moment  tensors  are  related  to  fault  plane  solutions. 
Body-wave  and  surface  wave  radiation  patterns  from  a 
source  represented  by  a  moment  tensor  are  discussed  by 
Kennett  (1988). 

Since  the  seismic  moment  tensor  Is  real  and  sym¬ 
metric,  a  principal  axis  transformation  can  be  found, 
diagonalizing  M  (Appendix  I).  The  diagonal  elements  are 
the  eigenvalues  m,  of  M.  Then,  the  scalar  seismic  moment 
can  be  determined  from  a  given  moment  tensor  by 

A/0  **  y  (  lmi  I  +  lm2 1  )  .  0°) 

where  m,  and  m2  are  the  largest  eigenvalues  (in  the  abso¬ 
lute  sense).  The  seismic  moment  can  equivalently  be 
estimated  by  the  relations  (Silver  and  Jordan,  1982): 


GENERAL  SEISMIC  POINT  SOURCES 

In  this  section,  it  is  assumed  that  the  seismic  source 
cannot  be  described  by  a  pure  double  couple  mechanism. 
The  moment  tensor  is  represented  as  sum  of  an  isotropic 
part,  which  is  a  scalar  times  the  identity  matrix,  and  a 
deviatoric  part. 

In  order  to  derive  a  general  formulation  of  the 
moment  tensor  decomposition,  let’s  consider  the  eigen¬ 
values  and  orthonormal  eigenvectors  of  the  moment  ten¬ 
sor.  Let  m,  be  the  eigenvalue  corresponding  to  the  ortho- 
normal  eigenvector  a,  =  (a,,  ,ait  .o,,)7".  Using  the  ortho¬ 
normality  of  the  eigenvectors  (Appendix  I,  (A1.5)),  we  can 
write  the  principal  axis  transformation  of  M  in  reverse 
order  as: 

M  -  [  a,  a2  a3  J  m 

°1»  a2l  a3l 
**  aly  a2y  a3 y 
aU  a2t  ®3i 

From  (21),  we  find  relations  between  components  of  the 
eigenvectors  and  moment  tensor  elements: 

Af„  -  mi  «?,  +  m2  a2l  +  m3  a32, 

Mn  -  m,  a,2,  +  m2  a 22,  +  m3  a32, 

Mu  -  m,  a  2,  +  m2  a 22,  +  m3  a32, 

M,,  -  "»i  «u  <»iy  +  "*2  a2l  a2,  +  m3  a3l  a3,  (22) 

Mu  “  ml  «li  ali  +  m2  ®2i  ®2j  +  m3  a3i  ®3* 

Mf,  —  m,  a a,,  +  m2  a2f  a2,  +  m3  a3y  a3l 

The  effect  of  the  eigenvalue  decomposition  (21)  is  that  a 

new  orthogonal  coordinate  system,  given  by  the  eigenvec¬ 
tors,  has  been  defined.  In  this  new  coordinate  system,  the 
source  excitation  is  completely  described  by  a  linear  com¬ 
bination  of  these  orthogonal  dipole  sources. 

m  in  (21)  is  the  diagonalized  moment  tensor.  The  ele¬ 
ments  of  m  are  the  eigenvalues  of  M.  We  now  define  the 


aV 


a3 


(21) 


m, 
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general  moment  tensor  decomposition  by  rewriting  m  as 


m 


tr(M)  0  0 

0  lr(M)  0 
0  0  <r(M) 


(23) 


+ 


m  i*  0  0 

0  m2  0 
0  0  m3 


3 


lr(M)  0 
0  lr(M) 
0  0 


0 

0 

tr(M) 


N 

+  £  m. 


i-i 


where  <r(M)  “  m  1  +  m2  +  m3  is  the  trace  of  the 
moment  tensor  and  m,  is  a  set  of  diagonal  matrices 
whose  sum  yields  the  second  term  in  (23).  The  purely 
deviatoric  eigenvalues  m‘  of  the  moment  tensor  are 


.  m  i  +  m2  +  m3  i 

m,  =  m, - - - -  m,  -  j  <r(M) 


(21) 


The  first  term  on  the  right  hand  side  (RHS)  of  (23) 
describes  the  isotropic  part  of  the  moment  tensor.  The 
eigenvalues  of  the  isotropic  part  of  the  moment  tensor  are 
important  for  quantifying  a  volume  change  in  the  source. 
The  second  term  describes  the  deviatoric  part  of  the 
moment  tensor  consisting  of  purely  deviatoric  eigenvalues, 
which  are  calculated  by  subtracting  1/3  tr  (M)  from  each 
eigenvalue  of  M.  This  deviatoric  part  of  the  moment  ten¬ 
sor  can  be  further  decomposed,  where  the  number  of 
terms  or  the  specific  form  of  the  decomposition  will  be 
discussed  in  the  next  sections.  Obviously,  a  multitude  of 
different  decompositions  are  possible.  In  Appendix  IN',  we 
give  some  numerical  examples  illustrating  several  methods 
of  moment  tensor  decomposition. 


Vector  Dipoles 

A  moment  tensor  can  be  decomposed  into  an  isotro¬ 
pic  part  and  three  vector  dipoles.  In  equation  (23)  let  N  = 
3  and 


m,"  0  0 

0  0  0 

0  0  0 

m,  - 

0  0  0 

0  0  0 

1X1  j  3 

0  m)  0 

0  0  0 

,  m, « 

0  0  0 

0  0 

, 

Applying  (21)  to  m,,  we  get  for  the  first  deviatoric  term 
(i=l)  in  the  decomposition 


Oil 

O  ll  0  ly 

“  li  0  )i 

lz  0  ly 

fl?y 

0  ly  0  1/ 

O 

liale 

aly“l/ 

an 

(20) 


where  we  identified  the  matrix  as  the  dyadic  a, a,  (Appen¬ 
dix  I).  The  dyadic  a,at  describes  a  dipole  in  the  direction 
of  the  eigenvector  a,.  By  applying  (21)  to  m2  and  m3  in 
(25),  we  get  similar  expressions  involving  a2&2  and  a3a3, 
describing  the  second  and  third  deviatoric  terms  in  the 
decomposition.  Finally,  equation  (21)  can  be  written  for 
the  decomposition  into  three  linear  vector  dipoles  along 
the  directions  of  the  eigenvectors  of  M  as 
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M  “  -^-(m,+m2-fm3)  I  (27) 

o 

+  m,  a,a,-|-m2  a2a2+m3  a3a3  , 

which  is  identical  to  (22)  and  equation  (4.55)  in  Ben- 
Menahem  and  Singh  (1081). 

Double  Couples 

Next,  we  decompose  a  moment  tensor  into  an  isotro¬ 
pic  part  and  three  double  couples.  For  the  deviatoric  part 
in  (23)  let  N  =  6  and 
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where  each  in,  is  equivalent  to  a  pure  double  couple 
source  (Appendix  III).  Notice  that  each  double  couple 
consists  of  two  linear  vector  dipoles  (c.f.  (25),  (26)  and 
(28)),  e.g.  (m,*/3)  (a, a,  —  a2a2)  for  in,.  Each  dipole  con¬ 
sists  of  two  forces  of  equal  strength  but  opposite  direction 
(c.f.  Figure  2).  Then,  the  double  couple  can  be  seen  this 
way:  The  first  couple  is  formed  by  one  force  of  each 
dipole,  one  force  pointing  in  the  positive  a,,  the  other  in 
the  negative  a2  direction.  The  corresponding  other  couple 
is  constructed  by  the  complementary  force  of  each  dipole, 
pointing  toward  the  negative  a,  and  positive  a2  direction. 
Using  (21)  with  (23)  and  (28),  we  get  the  result  that  a 
moment  tensor  can  be  decomposed  into  an  isotropic  part 
and  three  double  couples. 

M=  -i-(m  1+m2+m3)I+-i-(m1— m2)  (a,a,— a2a2)  (29) 

o  <5 

+  v(m2-m3)  (a2a2— a3a3)+^(m3— m,)  (a3a3-a,a,)  , 

o  6 

which  is  identical  to  equation  (4.57)  in  Ben-Menahem  and 
Singh  (1981). 


CLVD 

Alternatively,  a  moment  tensor  can  be  decomposed 
into  an  isotropic  part  and  three  compensated  linear  vector 
dipoles.  Adding  terms  like  in!  and  rn2  in  (28)  gives  a 
CLVD,  2a, a(  —  a2a2  —  a3a3.  This  CLVD  represents  a 
dipole  of  strength  2  in  the  direction  of  the  eigenvector  a,, 
and  two  dipoles  of  unit  strength  in  the  directions  of  the 
eigenvectors  a2  and  a3,  respectively.  The  decomposition 
can  then  be  expressed  as: 

M-  ~(m,+m2-l-m3)I+|-m1(2aIa1-a2a2-a3a3)  (30) 

+  jm^ajaj-a,  a, -a3a3)+y"»3(2a3a3-aia|-a2a2)> 


which  is  identical  to  equation  (4.56)  in  Ben-Menahem  and 
Singh  (1981). 


Major  and  Minor  Couple 

Next,  we  will  decompose  a  moment  tensor  into  an 
isotropic  component,  a  major  and  minor  double  couple. 
The  major  couple  seems  to  be  the  best  approximation  of  a 
general  seismic  source  by  a  double  couple  (Appendix  rV), 
since  the  directions  of  the  principal  axes  of  the  moment 
tensor  remain  unchanged.  The  major  double  couple  is  con¬ 
structed  in  the  following  way  (Kanamori  and  Given,  1981; 
Wallace,  1985):  The  eigenvector  of  the  smallest  eigenvalue 
(in  the  absolute  sense)  is  taken  as  the  null-axis.  Let's 
assume  that  |m3  |>|m2  |>|m,*  |  in  (23).  In  (23),  let 
N=2  and  use  the  deviatoric  condition  m,*+m2-t-m3  =  0 
to  obtain 
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Applying  (21)  to  m,t  we  get  the  first  deviatoric  term  in 
the  decomposition  which  corresponds  to  a  pure  double 
couple  termed  major  couple. 


Instead  of  the  major  double  couple,  a  best  double  couple 
can  be  constructed  similarly  by  replacing  m3  in  (32)  by 
the  average  of  the  largest  two  eigenvalues  (in  the  absolute 
sense,  Giardini,  1984).  Applying  (21)  to  m2  gives  the 
second  deviatoric  term  in  the  decomposition  which  also 
corresponds  to  a  pure  double  couple  termed  minor  couple. 


M*,,n  —  J  a,  a2  a3  J  m2 


(33) 


The  complete  decomposition  is  then: 


M= 


— (m  ,+m2+m3)I 


(34) 


+  m3  (a3a3-a2a2)  +ml*  (alal — *2*2) 


Double  Couple  -  CLVD 

Following  KnopofT  and  Randall  (1970)  and  Fitch  et 
al.  (1980),  we  can  decompose  a  moment  tensor  into  an  iso¬ 
tropic  part,  a  double  couple  and  a  compensated  linear  vec¬ 
tor  dipole.  Let's  assume  again  that  |m3  |>|m  2*  I  >  I  1*  I 
in  (23).  We  can  write  the  deviatoric  part  in  (23)  as  (N  — 
1) 


-F  0  0 

in,  -  m3  0  (F- 1)  0  ,  (35) 

0  0  1 

where  F  =  -  m|  /  m3  and  (F-l)  =*  m2  /  m3.  Note  that 
0<F  <0.5.  This  constraint  on  F  arises  from  the  deviatoric 
condition  m,*-l-m2-f-m3  =  0.  We  can  decompose  (35) 
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into  two  parts  representing  a  double  couple  and  a  CLVD 


0  0 

0 

[-1 

0 

0 

0  -1 

0 

-f  m3  F  0 

-1 
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O. 

1. 

1° 

0 
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where  we  assumed  that  the  same  principal  stresses 
produce  the  double  couple  as  well  as  the  CLVD  radiation. 
The  complete  decomposition  (21)  is  then: 

M“  j(m,+m2+m3)I  +  m3(l—2F)  (a^-a^)  (37) 

+  m3  F  (2a3a3-a2a2-alal) 


To  estimate  the  deviation  of  the  seismic  source  from 
the  model  of  a  pure  double  couple,  Dziewonski  tl  d. 
(1081)  used  the  parameter 


(38) 


where  m^in  is  the  smallest  eigenvalue  (in  the  absolute 
sense)  and  m^1)(  is  the  largest  (in  the  absolute  sense), 
given  by  (24).  From  (35),  we  see  that  c  =  F.  For  a  pure 
double  couple  source,  m„in  =  0  and  e  =  0;  for  a  pure 
CLVD,  (  =  0.5.  Alternatively,  t  can  be  expressed  in  per¬ 


centages  of  CLVD  (multiply  e  by  200.  The  percentage  of 
double  couple  is  (1— 2c)  *  100).  Dziewonski  and  Wood- 
house  (1983b,  see  also  Giardini,  1984)  investigated  the 
variation  of  c  versus  seismic  moment  and  earthquake  spa¬ 


tial  distribution  on  the  surface  of  the  earth. 


MOMENT  TENSOR  INVERSION 

There  are  various  methods  of  inversion  for  moment 
tensor  elements.  The  inversion  can  be  done  in  the  time  or 
frequency  domain.  Different  data  (e.g.  free  oscillations, 
surface-  and  body  waves;  different  seismogram  com¬ 
ponents)  can  be  used  separately  or  combined.  In  addition, 
certain  a  priori  constraints  such  as  tr  (M)  =  0,  or  M„  = 
Mf,  =  0  can  be  imposed  to  stabilize  the  inversion,  result¬ 
ing  in  a  decrease  in  number  of  resolved  moment  tensor 
elements.  In  this  Student’s  Guide,  we  briefly  outline  cer¬ 
tain  approaches  and  refer  to  the  original  papers  for 
further  reference. 

Gilbert  (1970)  introduced  the  seismic  moment  tensor 
for  calculating  the  excitation  of  normal  modes  (Saito, 
1967)  of  free  oscillations  of  the  earth.  Gilbert  (1973)  sug¬ 
gested  an  inversion  scheme  for  moment  tensor  elements  in 
the  frequency  domain.  Gilbert  and  Dziewonski  (1975) 
used  free  oscillation  data  for  their  moment  tensor  inver¬ 
sion.  Gilbert  and  Buland  (1976)  investigated  on  the  smal¬ 
lest  number  of  stations  necessary  for  a  successful  inversion 
(see  also  Stump  and  Johnson,  1977).  McCowan  (1976), 
Mendiguren  (1977),  Patton  and  Aki  (1979),  Patton  (1980), 
Romanowicz  (1981),  Kanamori  and  Given  (1981,  1982), 
Lay  el  af.  (1982),  Nakanishi  and  Kanamori  (1982,  1984), 
and  Scott  and  Kanamori  (1985)  used  long-period  surface 
waves  (typically  low  pass  filtered  at  135  sec).  Stump  and 
Johnson  (1977),  Strelitz  (1978,  1980),  Ward  (1980a,  b), 
Fitch  el  d.  (1980),  Langston  (1981),  Dziewonski  el  d 
(1981),  and  Dziewonski  and  Woodhouse  (1983a,  b),  used 
moment  tensor  inversion  for  body  wave  data  (typically 
low  pass  filtered  at  45  sec).  A  comparison  between 


moment  tensors  from  surface  waves  and  body  waves  was 
done  by  Fitch  el  d.  (1981).  Dziewonski  el  d  (1981)  sug¬ 
gested  an  iterative  inversion  method,  solving  for  the 
moment  tensor  elements  and  the  centroid  location  (Backus 
and  Mulcaby,  1076;  Backus,  1977a;  see  Dziewonski  and 
Woodhouse,  1083a  for  a  review).  The  reason  for  that 
approach  is  that  moment  tensor  elements  trade  off  with 
the  location  of  the  earthquake.  The  lateral  heterogeneity 
of  the  earth  was  considered  in  inversion  methods  by  Pat¬ 
ton  (1980),  Romanowicz  (1981),  Nakanishi  and  Kanamori 
(1982),  and  Dziewonski  el  of.  (1984c). 

The  moment  tensor  inversion  in  the  time  domain  can 
use  the  formulation  in  (5)  (e.g.  Gilbert,  1970;  McCowan, 
1976;  Stump  and  Johnson,  1977;  Strelitz,  1978;  Fitch 
el  d. ,  1980;  Ward,  1980b;  Langston,  1981).  If  the  source 
time  function  is  not  known  or  the  assumption  of  a  syn¬ 
chronous  source  is  dropped  (Sipkin,  1986),  the  frequency 
domain  approach  is  chosen  (e.g.  Gilbert,  1973;  Dziewonski 
and  Gilbert,  1974;  Gilbert  and  Dziewonski.  1975;  Gilbert 
and  Buland,  1976;  Mendiguren,  1977;  Stump  and  John¬ 
son,  1977;  Patton  and  Aki,  1979;  Patton,  1980;  Ward, 
1980a,  Kanamori  and  Given,  1981;  Romanowicz,  1981): 

dn(x,f)-Mki(f)Gnij(J)  (39) 

Both  approaches,  (5)  and  (39)  lead  to  linear  inversions  in 
the  time  or  frequency  domain,  respectively.  The  advan¬ 
tage  of  linear  inversions  is  that  a  large  number  of  fast 
computational  algorithms  are  available  (e.g.  Lawson  and 
Hanson,  1974;  Press  et  d  ,  1987).  We  can  write  either  (5) 
or  (39)  in  matrix  form: 

d  «  G  m  (40) 

In  the  time  domain,  the  vector  d  consists  of  n  sampled 
values  of  the  observed  ground  displacement  at  various 
arrival  times,  stations,  and  azimuths.  G  is  a  n  X  6 
matrix  containing  the  Green’s  functions  calculated  using 
an  appropriate  algorithm  and  earth  model,  and  m  is  a 
vector  containing  the  6  moment  tensor  elements  to  be 
determined  (Stump  and  Johnson,  1977).  In  the  frequency 
domain,  (40)  can  be  written  separately  for  each  frequency, 
d  consists  of  real  and  imaginary  parts  of  the  displacement 
spectra.  Weighting  can  be  introduced  which  actually 
smoothes  the  observed  spectra  subjectively  (Mendiguren, 
1977;  see  also  Ward,  1980b  for  weighting  of  body-wave 
data  in  the  time  domain).  In  the  same  way,  G  and  m 
contain  real  and  imaginary  parts,  m  contains  also  the 
transform  of  the  source  time  function  of  each  moment 
tensor  element.  If  constraints  are  applied  to  the  inversion, 
then  m  can  contain  a  smaller  number  of  moment  tensor 
elements.  In  such  a  case,  G  has  to  be  changed  accordingly. 
We  refer  to  Aki  and  Richards  (1980)  for  the  details  of 
solving  (40)  for  m  (Note  that  (40)  is  identical  to  their 
(12.83)). 

The  following  presents  an  outline  of  the  processing 
steps  in  a  moment  tensor  inversion.  The  first  step  is  the 
data  acquisition  and  the  preprocessing.  We  need  data 
with  good  signal  to  noise  ratio  that  are  unclipped  and 
that  have  a  good  coverage  of  the  focal  sphere  (Satake, 
1985).  Glitches  (non-seismic  high  amplitude  spikes  due  to 
non-linearity  of  instruments  e.g.  Dziewonski  ef  ol ,  1981) 
have  to  be  identified  and  possibly  removed.  Analog  data 
have  to  be  digitized.  The  effect  of  non-orthogonality  of 
the  analog  recorder  must  be  corrected.  The  digitized 
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record  has  to  be  interpolated  and  resampled  with  a  con¬ 
stant  sampling  rate.  At  this  point,  a  comparison  of  the 
sampled  waveform  with  the  original  one  can  help  to  iden¬ 
tify  digitization  errors.  The  horizontal  components  will 
be  rotated  into  radial  and  transverse  components.  Linear 
trends  have  to  be  identified  and  removed.  The  instrument 
effect  is  considered  next  (for  WWSSN  data  see  Hagiwara, 
1058;  for  SRO  data  see  McCowan  and  Lacoss,  1978;  for 
IDA  data  see  Agnew  et  ai,  1076).  We  can  use  either  one 
of  the  two  approaches;  i)  we  can  remove  the  instrument 
effect  from  the  observed  data  and  compare  with  theory  or 
ii)  we  can  apply  the  instrument  response  to  the  synthetic 
Green’s  functions  and  compare  with  observed  data.  The 
nominal  instrument  response  can  be  used  or  the  calibra¬ 
tion  of  the  instrument  can  be  checked  by  using  f.e.  the 
calibration  pulse  on  the  record.  In  addition,  the  polarity 
of  the  instruments  should  be  verified,  e.g.  from  records  of 
known  nuclear  explosions.  High  frequency  noise  in  the 
data  is  removed  by  low-pass  filtering.  Amplitudes  are 
corrected  for  geometrical  spreading  and  reflections  at  the 
free  surface  of  the  earth  (Bullen  and  Bolt,  1085).  For  sur¬ 
face  waves,  the  moving  window  analysis  (Landisman  et 
ai,  I960)  is  applied  in  order  to  determine  the  group  velo¬ 
city  dispersion.  From  this  analysis,  we  can  identify  the 
fundamental  mode  Rayleigh  and  Love  waves  which  can 
then  be  isolated. 

Second,  synthetic  Green’s  functions  are  calculated. 
Notice  that  the  Green’s  functions  are  dependent  on  the 
earth-model,  the  location  of  the  point  source  (centroid  of 
the  stress  glut,  or  epicenter  and  focal  depth),  and  the 
receiver  position. 

The  third  step  is  the  proper  inversion,  i.e.  the  solu¬ 
tion  of  (40)  (Aki  and  Richards,  1080).  Usually,  the  inver¬ 
sion  is  formulated  as  least  squares  problem  (Gilbert,  1073; 
Gilbert  and  Buland,  1976;  Mendiguren,  1977;  Stump  and 
Johnson,  1077).  However,  using  other  norms  can  have 
advantages  in  situations  where  less  sensitivity  to  gross 
errors  like  polarity  reversions  is  required  (Claerbout  and 
Muir,  1973;  Fitch  el  al.,  1080;  Patton,  1980). 

The  source  time  function  in  (5)  is  often  assumed  to 
be  a  step  function  (Gilbert,  1970,  1973;  McCowan,  1976; 
Stump  and  Johnson,  1077;  Patton  and  Aki,  1070;  Patton, 
1080;  Ward,  1080b;  Dziewonski  et  al.,  1081;  Kanamori 
and  Given,  1981).  Aiming  at  the  recovery  of  source  time 
functions,  Burdick  and  Mellman  (1076)  used  a  powerful 
iterative  waveform  inversion  method  based  on  optimizing 
the  cross-correlation  between  observed,  long-period  body- 
wave  trains  and  synthetics.  The  same  approach  was  used 
by  Wallace  et  al  ,  (1981)  in  order  to  invert  for  fault  plane 
solutions.  Other  methods  were  employed  by  Strelitz  (1980) 
and  Kikuchi  and  Kanamori  (1982)  for  large  earthquakes 
(see  also  Lundgren  et  al.  ,  1088).  Christensen  and  Ruff 
(1985)  reported  on  a  trade-off  between  source  time  func¬ 
tion  and  source  depth  for  shallow  events. 

If  the  focal  depth  is  not  known,  then  a  linear  inver¬ 
sion  can  be  done  for  each  depth  out  of  a  number  of  trial 
depths.  The  most  probable  depth  will  minimize  the  qua¬ 
dratic  error  between  observed  and  theoretical  waveforms 
(Mendiguren,  1077;  Patton  and  Aki,  1070;  Patton,  1080; 
Romanowicz,  1981).  The  influence  of  source  depth  on  the 
results  of  the  moment  tensor  inversion  was  investigated 
by  Sipkin  (1982;  Dziewonski  et  al. ,  1087b).  Differences  in 
source  depth  influence  the  relative  excitation  of  normal 


modes,  causing  systematic  errors  in  the  inversion. 

Systematic  errors  in  the  inversion  are  also  due  to 
deviations  of  the  earth-model  from  the  actual  properties 
of  the  earth,  affecting  the  synthetic  Green’s  functions. 
This  is  a  fundamental  problem  in  the  sense  that  we  are 
able  to  separate  the  source  effect  from  the  observed 
seismogram  only  to  a  limited  accuracy  (Mendiguren,  1977; 
Langston,  1081;  Silver  and  Jordan,  1982;  O'Connell  and 
Johnson,  1988).  A  major  problem  is  the  effect  of  lateral 
heterogeneity  of  the  earth  (Engdahl  and  Kanamori,  1080; 
Romanowicz,  1081;  Gomberg  and  Masters,  1088;  Sniedtr 
and  Romanowicz,  1088).  For  example,  a  relative  change 
of  0.5  %  due  to  lateral  heterogeneity  can  cause  a  misloca- 
tion  in  the  order  of  of  50  km  at  epicentral  distances  of 
about  00  degrees  (Dziewonski  and  Woodhouse,  1083b). 
Giardini  (1984)  and  Ekstrom  and  Dziewonski  (1985) 
reported  on  regional  shifts  in  centroid  positions  due  to 
lateral  heterogeneity.  In  the  inversion,  lateral  hetero¬ 
geneity  is  often  neglected,  i.e.  the  calculation  of  the 
Green’s  functions  is  usually  based  on  parallel  layers  of 
lateral  homogeneity  (Harkrider,  1064,  1070;  Langston  and 
Helmberger,  1075;  Harkrider,  1976).  Nakanishi  and 
Kanamori  (1082)  included  the  effect  of  lateral  hetero¬ 
geneity  into  the  moment  tensor  inversion.  Another 
approach  was  developed  for  earthquakes  within  a  small 
source  area:  a  calibration  event  is  declared  (mechanism 
known);  the  spectral  ratio  of  any  earthquake  in  that 
region  and  the  calibration  event  will  result  in  isolating  the 
difference  in  source  effects  -  the  influence  of  the  path  is 
eliminated  (Patton,  1080).  It  seems  that  the  errors  due  to 
lateral  heterogeneity  are  usually  large  enough  to  make  a 
statistical  significant  detection  of  an  isotropic  component 
of  the  moment  tensor  difficult  (Okal  and  Geller,  1979; 
Silver  and  Jordan,  1982;  Vasco  and  Johnson,  1988). 

Patton  and  Aki  (1070)  investigated  the  influence  of 
noise  on  the  inversion  of  long-period  surface  wave  data. 
They  found  that  additive  noise  such  as  background 
recording  noise  does  not  severely  affect  the  results  of  a 
linear  inversion.  However,  multiplicative  noise  (signal 
generated  noise)  caused  by  focusing,  defocusing,  mul- 
tipathing,  higher  mode  or  body  wave  interference,  and 
scattering  distorts  the  inversion  results  significantly 
(overestimation  or  underestimation  of  moment  tensor  ele¬ 
ments,  deviation  from  the  source  mechanism;  Patton, 
1080;  Ward,  1980b).  Finally,  body  waves  of  events  with 
moments  larger  than  1027  dyne-cm  are  severely  affected  by 
finiteness  of  the  source  and  directivity.  If  not  corrected 
for,  an  inversion  can  lead  to  severe  errors  in  the  moment 
tensor  elements  (Dziewonski  et  ai,  1981;  Kanamori  and 
Given,  1081;  Patton  and  Aki,  1070;  Lay  et  ai,  1982; 
Giardini,  1084). 

The  inversion  has  only  a  limited  resolution  of 
moment  tensor  elements  for  certain  data.  If  we  have  spec¬ 
tra  of  fundamental  mode  Rayleigh  waves  only,  the  con¬ 
straint  that  the  trace  of  the  moment  tensor  vanishes  (no 
volume  change)  must  be  applied  (Mendiguren,  1077;  Pat¬ 
ton  and  Aki,  1979).  This  constraint  is  linear.  Another 
constraint  which  is  often  applied  in  addition  is  that  one 
eigenvalue  vanishes  (approximating  the  source  by  a  dou¬ 
ble  couple).  This  constraint  is  not  linear  (Strelitz,  1978: 
Ward,  1080b).  In  such  a  case,  the  inversion  is  iterative, 
using  a  linearized  version  of  the  constraints  (Ward, 
1080b).  For  earthquakes  at  shallow  depths  (less  than 
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about  30  km),  the  moment  tensor  elements  Mtt  and 
corresponding  to  vertical  dip  slip  faulting  are  not  well 
constrained  from  long-period  surface  wave  data  since  the 
related  excitation  functions  assume  very  small  values  near 
the  surface  of  the  earth  (Fitch  et  ai,  1981;  Dziewonski  et 
a/.,  1981;  Kanamori  and  Given,  1981,  1982;  Dziewonski 
and  Woodhouse,  1983a).  In  order  to  overcome  this  prob¬ 
lem,  additional  independent  data,  such  as  fault  strike 
(observed  surface  breakage)  can  be  introduced  into  the 
inversion.  Another  approach  is  to  constrain  these 
moment  tensor  elements  to  be  zero.  Thus,  possible  fault 
mechanisms  are  restricted  to  vertical  strike  slip  or  45 
degree  dip  slip  (Kanamori  and  Given,  1981,  1982). 

In  Appendix  V,  we  relate  the  Green’s  functions  in  the 
formulation  of  Herrmann  and  Wang  (1985)  to  a  simple 
moment  tensor  inversion  scheme.  This  inversion  example 
is  aimed  at  testing  computer  programs. 

CONCLUSION 

A  seismic  moment  tensor  describes  the  equivalent 
forces  of  a  seismic  point  source.  The  eigenvectors  are  the 
principal  axes  of  the  seismic  moment  tensor.  For  pure 
double  couple  sources,  the  principal  axis  corresponding  to 
the  negative  eigenvalue  is  the  pressure  axis,  the  principal 
axis  corresponding  to  the  positive  eigenvalue  is  the  tension 
axis,  and  the  principal  axis  corresponding  to  the  eigen¬ 
value  zero  gives  the  null  axis.  The  pressure,  tension,  and 
null  axes  can  be  displayed  in  the  familiar  focal  mechanism 
plot  (fault  plane  solution).  For  general  seismic  sources,  we 
can  decompose  the  seismic  moment  tensor.  First,  we  can 
separate  out  the  isotropic  component  which  describes  the 
volume  change  in  the  source.  The  leftover  part  of  the 
moment  tensor  has,  in  general,  three  nonvanishing 
eigenvalues.  This  deviatoric  part  of  the  moment  tensor 
can  be  decomposed  into  a  number  of  simple  combinations 
of  equivalent  forces.  Obviously,  there  is  no  unique 
moment  tensor  decomposition,  i.e.  unique  model  of 
equivalent  forces.  We  outlined  methods  of  determining 
moment  tensor  elements  from  observations,  indicating 
that  recording  noise  as  well  as  systematic  errors  due  to  an 
insufficient  knowledge  of  the  Green’s  functions  can  intro¬ 
duce  errors  into  the  moment  tensor  elements.  This  sug¬ 
gests  caution  when  apparent  non-double  couple  sources 
result  from  the  inversion. 

Randall  and  Knopoff  (1970),  Gilbert  and  Dziewonski 
(1975),  Dziewonski  et  at.  (1981),  Kanamori  and  Given 
(1981,  1982),  Dziewonski  and  Woodhouse  (1983b),  Giar- 
dini  (1984),  and  Scott  and  Kanamori  (1985)  reported  that 
some  seismic  sources  cannot  be  described  by  a  pure  double 
couple.  One  explanation  is  that  some  fault  planes  show  a 
complex  geometry  (Dziewonski  and  Woodhouse,  1983b). 
Another  explanation  can  be  that  some  sources  deviate 
from  the  model  of  a  sudden  shear  dislocation;  they  can  be 
due  to  a  rapidly  propagating  phase  transition  (Knopoff 
and  Randall,  1970;  Dziewonski  and  Woodhouse,  1983b). 
However,  the  simple  inversion  experiment  in  Appendix  V 
pointed  out  that  the  deviation  from  a  pure  double  couple 
can  also  be  due  to  the  presence  of  noise  in  the  data 
(Stump  and  Johnson,  1977;  Patton  and  Aki,  1979;  Pat¬ 
ton,  1980;  Ward,  1980b;  Wallace,  1985;  O'Connell  and 
Johnson,  1988). 
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APPENDIX  I 

In  the  following,  we  give  some  mathematical 
definitions  of  tensors,  the  eigenvalue  problem  and  dyadics 
following  Arfken  (1985). 

Let  M  be  a  moment  tensor  of  second  rank  (order). 
Then,  M  is  represented  as  a  3X3  matrix  in  a  given  refer¬ 
ence  frame.  Let  afk  be  the  cosine  of  the  angle  between  the 
p  axis  of  another  coordinate  system  and  the  k  axis.  Then 
the  components  of  M,  Mk] ,  transform  into  the  new  refer¬ 
ence  frame  by  the  relation 

“  akP  a„  •  (A  1.1) 

where  we  need  to  sum  over  repeated  indices  (summation 
convention). 

Given  a  moment  tensor  M,  let’s  assume  that  there  is 
a  vector  a  and  a  scalar  m  such  that 

M  a  =  m  a  (A1.2) 

a  is  called  eigenvector  of  M  and  m  is  the  corresponding 
eigenvalue.  For  calculating  the  eigenvalues  and  eigenvec¬ 
tors  of  a  given  moment  tensor  (solving  the  eigenvalue 
problem),  we  transform  (A1.2) 

(M-ml)a-O  ,  (A1.3) 

where  I  is  the  identity  matrix.  (A1.3)  is  a  system  of  3 
simultaneous  homogeneous  linear  equations  in  at.  Non¬ 
trivial  solutions  are  found  by  solving  the  secular  equation 
(characteristic  polynomial) 

det(M  —  ml)  =  0  ,  (A1.4) 

where  "det"  means  the  determinant.  (A1.4)  is  a  polyno¬ 
mial  of  third  degree.  It  has  three  real  roots,  i.e.  eigen¬ 
values,  since  the  moment  tensor  is  real  and  symmetric 
(Faddeeva,  1959).  Substituting  each  eigenvalue  m,  into 
(A1.3)  gives  the  corresponding  eigenvector  a,.  The  eigen¬ 
vectors  are  orthogonal.  Multiplying  each  eigenvector  by 
its  inverse  norm,  we  get  the  orthonormal  eigenvectors, 
renaming  them  as  a,-: 

»<*;“*«  •  (Al  .5) 

Knowing  the  eigenvectors,  we  can  diagonalize  M  (princi¬ 
pal  axis  transformation).  Let  A  be  the  matrix  whose 
columns  are  the  orthonormal  eigenvectors  of  M.  From 
(A1.5),  we  see  that  A  is  orthogonal  :  A T  =  A-1.  Then. 
ArMA-*m,  where  m  is  diagonal,  consisting  of  the 
eigenvalues  of  M. 

We  represent  a  dyadic  by  writing  two  vectors  a  and 
b  together  as  ab  (see  Appendix  A  in  Ben-Menahcm  and 
Singh,  1981).  These  two  vectors  forming  the  dyadic  are 
not  operating  on  each  other,  but  define  a  3X3  matrix.  Let 
I,  j,  and  k  be  unit  vectors  along  a  right  handed  Cartesian 
coordinate  system.  The  dyadic  ab  is  defined  as 
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ab  -  (a,i  +  a, j  +  axk)  (4xi  -f  4(  j  +  4,k) 

—  iia,  b,  +  ij ax  bf  +  ika,  4, 

+  jia,6I  +  jja,4,  +  jka,4l 
+  kia2ft,  +  kjax4f  +kka,6I  (A1.6) 

a 1 4t  «*  4,  o,  4, 

“  ajr  °jr  ^jr  fl|r 

ai  4j  a,  bf  a,  b. 

For  a  =  b,  we  get  (26).  The  multiplication  of  a  vector  c 
from  the  left  is 

c  ij  =  [(ict  +  j cf  +  kcj  i]j  =  cj  (A1.7) 

If  the  dyadic  is  symmetric,  the  multiplication  of  any  vec¬ 
tor  with  the  dyadic  is  commutative,  i.e.  ab  =  ba.  In 
general,  we  can  understand  a  dyadic  as  a  tensor  of  second 
rank.  By  a  proper  choice  of  the  coordinate  system,  a  sym¬ 
metric  dyadic  can  always  be  transformed  into  diagonal 
form  (principal  axis  transformation).  As  an  example,  we 
can  rewrite  (10)  using  dyadics  (Gilbert,  1973): 

ui/  +  Mi  *=  tt  —  pp  (A1.8) 

*  0.5  l(t+p)(t-p)  +  (t-p)(t+p)]  . 


APPENDIX  II 

The  PDE  monthly  listings  (NEIS)  routinely  publish 
centroid  moment  tensor  solutions  in  the  notation  of  the 
normal  mode  theory  following  Dziewonski  et  al  (1981). 
For  reference,  the  spherical  moment  tensor  elements,  /,  , 
in  the  notation  of  Gilbert  and  Dziewonski  (1975), 
Dziewonski  et  al.  (1981),  and  Dziewonski  and  Woodhouse 
(1983a)  are  compared  to  the  moment  tensor  elements  as 
given  in  (18)  following  Aki  and  Richards  (1980). 


/. 
f  2 

!  3 

f* 

J  5 
/ft 


M„  ■ 

“  Af„ 

^ee 

=  M„ 

=  Mn 

A're 

K* 

Me* 

1 

1 

£ 

(A2.1) 


where  (r,9,<4)  are  the  geographical  coordinates  at  the 
source.  0  is  the  colatitude  (0  =  0  at  the  north  pole)  and 
<j>  is  the  longitude  of  the  point  source.  The  sign  of  the  off- 
diagonal  moment  tensor  elements  depend  on  the  orienta¬ 
tion  of  the  coordinate  system  (Fitch  el  al.,  1981).  But  the 
eigenvalues  and  the  eigenvectors  of  the  moment  tensor  in 
the  formulation  of  (18)  or  (A2.1)  are  identical,  which  can 
be  shown  by  comparing  the  solutions  to  the  secular  equa¬ 
tion  (Appendix  1).  This  result  is  expected  since  physical 
laws  should  not  depend  on  the  choice  of  the  reference 
frame.  The  slip  vector  u  and  fault  normal  v  are 
(Dziewonski  and  Woodhouse,  1983a) 


vi  -■  ti  (—  cos  X  cos  ♦  —  cos  5  sin  X  sin  $  )  ee 


+  u  (  cos  X  sin  ♦  —  cos  6  sin  X  cos  *  )  e,  (A2.2) 

+  #  sin  X  sin  4  er  , 


v  —  sin  S  sin  4>  ee  +  sin  S  cos  -f  cos  S  er  (A2.3) 

These  two  equations  are  identical  to  (4.122)  in  Ben- 
Menahem  and  Singh  (1981).  The  differences  in  sign  com¬ 
pared  to  (15)  and  (18)  can  be  fully  explained  by  noting 
that  er  =  -  e, ,  =  e( ,  and  ee  =  -  ex ;  in  other  words, 

ee,  e^,  and  er  are  unit  vectors  towards  south,  east,  and 
up,  respectively  (defining  a  right  handed  system). 

APPENDIX  HI 

In  order  to  gain  some  experience  in  the  relationships 
between  a  moment  tensor  and  a  fault  plane  solution,  three 
simple  focal  mechanisms  are  discussed  in  detail.  These 
will  be  vertical  strike  slip,  45  degree  dip  slip,  and  vertical 
dip  slip  faults.  These  three  fault  plane  solutions  form  a 
complete  set  :  The  seismic  radiation  from  a  dislocation  on 
a  plane  dipping  an  arbitrary  angle  (but  striking  north- 
south)  can  be  expressed  as  a  linear  combination  of  these 
three  solutions  (Burridge  et  al,  1964;  Ben-Menahem  and 
Singh,  1968). 


Vertical  strike  slip  fault 


The  following  focal  mechanism  is  assumed:  (strike)  ^ 


—  0°,  (dip)  6  =  90°,  and  (slip)  X  =  0°.  From  (15)  and 
(16),  the  slip  vector  on  the  fault  plane  is  u  =  (1,0,0)  and 
the  vector  normal  to  the  fault  plane  is  v  =  (0,1,0).  The 
moment  tensor  can  be  determined  from  (18). 


M  - 


0  M0  0 
M0  0  0 
0  0  0 


(A3.1) 


The  eigenvalues  and  eigenvectors  of  this  tensor 
(Example  4.6.1  in  Arfken,  1985,  see  also  Appendix  I)  are 
shown  in  Table  A.l  (The  components  of  the  eigenvectors 
are  north,  east,  and  down). 


TABLE  A.l 


EIGENVALUE 

EIGENVECTOR 

0 

(  0.0000,  0.0000,-1.0000) 

M0 

(  0.7071,  0.7071,  0.0000) 

■Ma 

(-0.7071,  0.7071,  0.0000) 

The  eigenvector  b  corresponding  to  the  eigenvalue 
zero  is  the  null-axis,  the  eigenvector  t  corresponding  to 
the  positive  eigenvalue  gives  the  tension  axis,  T,  and  the 
eigenvector  p  corresponding  to  the  negative  eigenvalue 
gives  the  pressure  axis,  P,  of  a  focal  mechanism. 

The  focal  mechanism  is  obtained  by  using  (7)-(14) 
(Herrmann,  1975).  For  the  trend  and  plunge  (in  degrees) 
of  the  X-,  Y-,  null-,  T-,  and  P-axes,  we  get  (90,  0),  (180, 
0),  (270,  90),  (45,  0),  and  (135,  0),  respectively.  The  trend 
of  both  the  P  and  T  axes  can  be  shifted  by  180r  (Figure 
A. la);  i.e.  the  P-axis  can  also  be  described  by  (315*,  0") 
and  the  T-axis  by  (225*,  0° ).  This  ambiguity  can  be  fol¬ 
lowed  through  to  the  moment  tensor:  The  sign  of  an 
eigenvector  is  not  constrained  by  the  solution  of  the  eigen¬ 
value  problem  (Arfken,  1985).  However,  any  choice  of 
sign  leads  to  the  same  focal  mechanism. 


and 
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Fig.  A.l.  Focal  mechanisms  of  a  vertical  strike  slip  fault  (strike  ==  0°,  dip 
=  90° ,  slip  =  0° ),  (a),  a  45  degree  dip  slip  fault  (strike  =  0° ,  dip 
=  45°,  slip  =  90° ),  (b),  and  a  vertical  dip  slip  fault  (strike  —  0°, 
dip  =  90°,  slip  =  90°),  (c)  (Appendix  III). 


45  degree  dip  slip  fault 


The  following  focal  mechanism  is  assumed:  (strike)  4> 
—  0*,  (dip)  6  =  45°,  and  (slip)  X  —  90°.  From  (15)  and 
(16),  u  =  (0,-0.7071,-0.7071)  and  u  —  (0,0.7071,-0.7071). 
The  moment  tensor  is  calculated  from  (18). 


M  - 


0  0  0 
0  — Af0  0 
0  0  M0 


(A3.2) 


The  corresponding  eigenvalues  and  eigenvectors  are 
shown  in  Table  A.2. 


TABLE  A.2 


eigenvalue 

eigenvector 

0 

(-1, 0, 0) 

M0 

(  0,  0,-1) 

- zMn _ 

_ (  0.  1.  0) 

The  fault  plane  solution  is  obtained  from  (7)-(14) 
(Herrmann,  1975).  For  the  trend  and  plunge  (in  degrees) 
of  the  X-,  Y-,  null-,  T-,  and  P-axes,  we  get  (90,  45),  (270, 
45),  (360,  0),  (180,  00),  and  (270,  0),  respectively.  The 
trend  of  the  P  and  null  axes  can  be  shifted  by  180°  (Fig¬ 
ure  A.lb)  to  (90°  ,  0°  )  and  (180°  ,  0° ),  respectively. 


Vertical  dip  slip  fault 


The  following  focal  mechanism  is  assumed:  (strike)  <f> 
==  0°,  (dip)  5  =  90°,  and  (slip)  X  =  90°.  From  (15)  and 
(16),  u  =  (0,0,-l)  and  u  =  (0,1,0).  The  moment  tensor  is 
calculated  from  (18). 


M  - 


0  0  0 
0  0  ~M0 

o  —m0  o 


(A3.3) 


The  corresponding  eigenvalues  and  eigenvectors  are 
shown  in  Table  A.3. 
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TABLE  A.3 


EIGENVALUE 

EIGENVECTOR 

0 

(-1.0000,  0.0000,  0.0000) 

A/0 

(  0.0000,  0.7071,-0.7071) 

_ zMn _ 

(  0.0000.  0.7071,  0.7071) 

The  fault  plane  solution  is  obtained  from  (7)-(14) 
(Herrmann,  1075).  For  the  trend  and  plunge  (in  degrees) 
of  the  X-,  Y-,  null-,  T-,  and  P-axes,  we  get  (0.  00).  (00,  0). 
(180,  0),  (270,  45),  and  (00,  45),  respectively.  The  trend 
of  the  null  axis  can  be  shifted  by  180°  (Figure  A.lc)  to 
(3G0° ,  0° ). 


APPENDIX  IV 


In  the  following,  examples  of  the  five  methods  of 
moment  tensor  decomposition  are  presented. 


In  order  to  construct  a  moment  tensor  that  does  not 
lead  to  a  simple  double  couple  mechanism,  let 


1 

0 

0 

M, 

-  1 

0 

1 

0 

.0 

0 

1 

. 

0 

I 

0 

m2 

=  6 

1 

0 

0 

0 

0 

0 

0 

0 

0 

m3 

=  3 

0 

-1 

0 

0 

0 

1 

0 

0 

0 

-  1 

0 

0 

— 

.0 

-1 

0 

(A4.1) 


(A4.2) 


(A4.3) 


(A4.-I) 


The  first  moment  tensor  represents  an  explosion,  the 
others  arc  the  familiar  ones  from  Appendix  III,  represent¬ 
ing  a  vertical  strike-slip,  a  45  degree  dip-slip,  and  a  verti¬ 
cal  dip-slip  fault,  respectively.  All  four  moment  tensors 
are  superimposed  in  order  to  describe  a  complex  source 
that  is  dominated  by  a  vertical  strike  slip  mechanism. 
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Fig.  A.2.  Focal  mechanisms  of  the  double  couples  from  the  moment  tensor 
decomposition  (Appendix  IV).  (a)  major  couple  of  the  moment  ten¬ 
sor  in  (A4.5),  elementary  moment  tensor  EMT3  in  (A4.6),  and 
second  term  on  the  RHS  of  (A4.9)  (strike  =  355",  dip  =  80",  slip 
—  16"),  (b)  elementary  moment  tensor  EMT2  in  (A4.6)  (strike  = 
125",  dip  =  63°,  slip  =  —05"),  (c)  elementary  moment  tensor 
EMT4  in  (A4.6)  (strike  =  199",  dip  =  44",  slip  =  63"). 


The  result  is 


M 


1  6  0 
6  -2  -1 
0-14 


(A4.5) 


Table  A. 4  shows  the  eigenvalues  of  (A4.5)  and  the 
corresponding  eigenvectors,  which  are  the  principal  axes  of 


M. 


TABLE  A.4 


EIGENVALUE 

EIGENVECTOR 

3.8523 

5.8904 

-6.7427 

(-0.2938,  -0.1307,  -0.9456) 
(  0.7352,  0.5002,  -0.3170) 

(  0.6109.  -0.7883,  -0.07341 

The  sum  of  the  eigenvalues  is  equal  to  3,  which  is  the 
expected  value  for  the  sum  of  the  eigenvalues  of  (A4.1), 
describing  an  explosion. 

In  order  to  calculate  the  deviatoric  part  of  the  given 
moment  tensor,  the  isotropic  part  is  removed  by  subtract¬ 
ing  one  third  of  the  trace  of  (A4.5)  from  each  diagonal  ele¬ 
ment.  The  solution  to  the  corresponding  eigenvalue  prob¬ 
lem  leads  to  the  same  eigenvectors  as  above.  This 
indicates  that  the  principal  axes  of  the  complete  moment 
tensor  are  the  same  as  the  principal  axes  of  the 
corresponding  deviatoric  tensor.  The  deviatoric  eigen¬ 
values  are  2.8523,  4.8904,  and  -7.7427  in  the  order  of 
Table  A.4  (see  (24)).  From  (38),  t  —0.37,  i.e.  the  given 
moment  tensor  has  a  double  couple  component  of  26  % 
and  a  CLVD  component  of  74  %. 

F»r  the  determination  of  the  major  couple  from  (32), 
we  identify  the  eigenvector  (0.6109,  -0.7883,  -0.0734) 
corresponding  to  the  deviatoric  eigenvalue  of  -7.7  as  the 
P-axis,  the  eigenvector  (0.7352,  0.5992,  -0.3170) 

corresponding  to  the  eigenvalue  of  4.9  as  the  T-axis,  and 


the  eigenvector  (-0.2938,  -0.1397,  -0.9456)  with  eigenvalue 
2.9  as  the  null-axis.  The  fault  plane  solution  of  the  major 
double  couple  gives  for  the  X-,  Y-,  null-,  T-,  and  P-axes 
(in  degrees):  (172,  16),  (265,  10),  (25,  71),  (219,  18),  and 
(128,  4),  respectively  (Figure  A.2a).  The  major  double 
couple  gives  a  good  estimate  of  the  major  contribution  to 
the  faulting  which  is  predominar  ly  strike  slip  (compare 
Figures  A. la  and  A. 2a). 

Next,  the  moment  tensor  in  (A4.5)  is  decomposed 
into  an  isotropic  part  and  three  double  couples  following 
(29)  which  is  evaluated  by  using  (26)  together  with  the 
data  in  Table  A.4.  The  numbering  of  the  eigenvalues  and 
eigenvectors  in  (29)  follows  the  columns  of  Table  A.4,  but 
that  is  not  relevant  to  the  solution.  The  calculation  gives 


1  0  0 

0  4542  0  3995  -0.5109 

M  =  1 

0  1  0 

+  0.6794 

0.3995  0  3396  -0.3220 

0  0  1 

-0  5109  -0  3220  -0.7936 

+  4.2110 


+  3.5316 


0  1673 
09221 
-0.1882 


0.9221 
-0  2623 
-02477 


-0.1882 

-0.2477 

0.0951 


-02869  0.5226  0.3227 
0  5226  -0  6019  0.0743 
0.3227  0.0743  0.8887 


(A4.6) 


This  equation  is  identical  to  (A4.5).  The  first  ele¬ 
mentary  moment  tensor  (EMT1)  on  the  RHS  of  (A4.6) 
describes  the  explosion  (isotropic  component  of  (A4.5)) 
and  is  identical  to  (A4.1).  The  last  three  elementary 
moment  tensors  on  the  RHS  (EMT2,  EMT3,  EMT4, 
respectively)  represent  pure  double  couple  sources  since 
the  eigenvalues  of  each  tensor  is  0  and  ±  1.  The  three  ele¬ 
mentary  moment  tensors  have  identical  eigenvectors 
which  are  the  same  vectors  as  shown  in  Table  A.4. 
However,  the  correlation  between  eigenvector  and  eigen¬ 
value  (i.e.  null-,  P-,  and  T-axes)  varies.  Note  that  replac- 
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ing  Mtj  by  —  A/t;-  switches  the  sign  of  the  eigenvalues 
(leaving  the  eigenvectors  untouched),  which  is  identical  to 
interchanging  the  P-  and  T-axes. 

From  the  eigenvalues  and  the  eigenvectors,  the  fault 
plane  solution  for  each  elementary  moment  tensor  is 
determined  and  shown  in  Table  A.5. 


TABLE  A.5 


AXIS 

EMT2 

TRD 

(deg.) 

EMT2 

PLG 

(deg.) 

EMT3 

TRD 

(deg.) 

EMT3 

PLG 

(deg.) 

EMT4 

TRD 

(deg.) 

EMT4 

PLG 

(deg.) 

X 

36 

26 

172 

16 

324 

38 

Y 

226 

63 

265 

10 

109 

46 

NULL 

128 

4 

25 

71 

219 

18 

T 

219 

18 

219 

18 

25 

71 

P 

25 

71 

128 

4 

128 

4 

The  focal  mechanisms  corresponding  to  EMT2  - 
EMT4  are  shown  in  Figures  A.2b,  A.2a,  and  A.2c,  respec¬ 
tively.  Note  that  the  positions  of  the  axes  remain  fixed  in 
these  figures,  where  only  the  correlation  to  the  eigenvalues 
changes.  The  fault  plane  solution  representing  the  third 
elementary  moment  tensor  EMT3  in  (A4.6)  is  identical  to 
the  fault  plane  solution  of  the  major  couple  (Figure  A.2a). 
Notice  that  this  solution  has  also  the  largest  coefficient  in 
(A4.6).  This  solution  is  an  approximation  to  the  major 
contributor  of  the  moment  tensor  (Figure  A.la  and 
(A4.2)).  However,  the  other  fault  plane  solutions  (Figure 
A.2b  and  A.2c)  do  not  show  similarities  to  the  input  fault 
mechanisms  (Figure  A.lb  and  A.lc). 

The  seismic  moments  of  the  elementary  moment  ten¬ 
sors  are  given  by  the  coefficients  in  (A4.6).  The  sum  of  the 
seismic  moments  of  the  elementary  moment  tensors  is  1.4 
times  larger  than  the  seismic  moment  of  the  composite 
moment  tensor  in  (A4.5). 

Next,  the  moment  tensor  in  equation  (A4.5)  is 
decomposed  into  an  isotropic  part  and  three  vector  dipoles 
following  (27)  which  is  evaluated  by  using  (20)  together 
with  Table  A.4. 


1  0  0 

0  0863  0  0410  0.2779 

M  -  1 

0  1-0 

+  2.8523 

0  0410  0.0195  0.1321 

.0  0  1 

.0.2779  0.1321  0.8941 

+  4  8904 

-  7.7427 


0.5405 
0  4405 
-0  2330 
0.3732 
-0  4818 
-0  0448 


04405 

03591 

-0.1899 

-0.4816 

0.6214 

0.0578 


-0.2330 

-0.1899 

01005 

-0.0448 

0.0578 

00054 


(A4.7) 


This  equation  is  identical  to  (A4.5).  In  the  notation 
used  above,  each  of  the  elementary  moment  tensors 
EMT2,  EMT3,  and  EMT4  have  two  eigenvalues  equal  to 
xero,  the  third  one  equals  one.  EMT2  is  represented  by  the 
eigenvector  (0.2938,  0.1307,  0.0456),  EMT3  by  (-0.7352, 
-0.5092,  0.3170),  and  EMT4  by  (0.6100,  -0.7883,  -0.0734). 
These  vector  dipoles  are  mutually  orthonormal.  Notice 
that  these  vector  dipoles  are  identical  to  the  eigenvectors 
of  equation  (A4.5),  which  are  the  principal  axes  of  the  ten¬ 
sor  (Table  A.4).  EMT2  represents  the  null-,  EMT3  the 
tension-,  and  EMT4  the  pressure  axis.  The  seismic 


moments  of  the  elementary  moment  tensors  are  given  by 
the  coefficients  in  (A4.7)  which  are  identical  to  the  devia- 
toric  eigenvalues  of  (A4.5).  This  exercise  demonstrated 
that  vector  dipoles  are  related  to  the  eigenvectors  scaled 
by  the  corresponding  eigenvalue  of  a  given  moment  ten¬ 
sor,  which  makes  an  evaluation  of  (A4.7)  obsolete. 

Alternatively,  the  moment  tensor  in  equation  (A4.5) 
can  be  decomposed  into  an  isotropic  part  and  three  com¬ 
pensated  linear  vector  dipoles  using  (30). 


1  0  0 

-07411  0  1231  08336 

M  =  1 

0  1  0 

+  1.2841 

0.1231  -0.9415  0.3963 

0  0  1 
f 

0  8336  0  3963  1.6823 

■a 

0.6215  1.3216  -0  6991 


+  1.9635  1.3216  0  0773  -0  5697 

-0.6991  -0.5697  -0.6985 
0.1196  -1  4447  -0  1345 
-  2.2476  —1.4447  0.8642  0.1734 

-0.1345  0.1734  -0  9838 


(A4.8) 


This  equation  is  identical  to  (A4.5).  The  seismic 
moments  of  the  elementary  moment  tensors  are  given  by 
the  product  of  the  respective  coefficient  and  V3  .  The 
eigenvalues  and  eigenvectors  for  (A4.8)  are  shown  in 
Table  A.6,  using  the  same  notation  as  above.  Note  that 
the  eigenvectors  are  identical  to  those  in  Table  A.4. 


TABLE  A.6 


EMT 

EIGENVALUE 

EIGENVECTOR 

2 

-1 

(  0.6109,-0.7883,-0.0734) 

-1 

(  0.7352,  0.5992,-0.3170) 

2 

(-0.2938,-0  1397.-0.9456) 

3 

-1 

(-0.2938,-0.1307,-0.0456) 

-1 

(  0.6109,-0.7883,-0.0734) 

2 

(  0.7352,  0.5992.-0.3170) 

4 

-1 

(  0.7352,  0.5992,-0.3170) 

-1 

(-0.2938,-0.1307,-0.0456) 

2 

(  0.6109,-0.7883.-0.0734) 

Next,  the  moment  tensor  in  (A4.5)  is  decomposed 
into  an  isotropic  part,  a  double  couple  and  CL\T>  follow¬ 
ing  (37),  where  e  =  F  =  0.3684. 


I  0  0 

0.1673  09221  -0.1882 

M  =  1 

0  1  0 

+  2.0379 

0.9221  -0.2623  -0.2477 

0  0  1, 

-0.1882  -0.2477  0  0951 

-  2.8523 


0  1196  -1.4447  -0.1345 

-1.4447  0.8642  0.1734 

,-0.1345  0  1734  -0.9838 


(A4.9) 


This  equation  is  identical  to  (A4.5).  Notice  that  the 
second  term  on  the  RHS  corresponds  to  EMT3  in  (A4.6) 
and  to  the  major  double  couple.  These  three  tensors  all 
have  the  same  fault  plane  solution  (Figure  A.2a).  The 
third  term  in  (A4.9)  corresponds  to  EMT4  in  (A4.8), 
representing  a  CLVD  (see  Table  A.6). 

As  final  remark,  let's  consider  the  decomposition 
equations  (27),  (29),  and  (30)  for  a  simple  double  couple 
source  (<  =  0  ),  e.g.  let  m j  =  -  m2  =  1,  and  m3  =  0. 
Then,  M  =  -  a2a3  for  all  three  equations.  That  is, 

we  get  one  pure  double  couple  out  of  the  decomposition. 
For  a  CLVD  (e  =  0.5),  let’s  assume  that  nq  =  m2  =  -1, 
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and  m3  =  2.  Then  all  three  formulas  give  M  =  2  a3a3  - 
a,aj  -  a2a2,  representing  one  CLVD. 

APPENDIX  V 

In  this  section,  we  relate  the  Green's  functions  in  the 
formulation  of  Herrmann  and  Wang  ^1985)  to  a  moment 
tensor  inversion  scheme.  Following  the  theory  given  by 
Herrmann  and  Wang  (1985),  the  Fourier  transformed  dis¬ 
placements  at  the  free  surface  at  the  distance  r  from  the 
origin  due  to  an  arbitrarily  oriented  double  couple 
without  moment  is 

</2(r,z=0,w)  =  ZSS  A  i  +  ZDS  A2  +  ZDD  A3 

dr(r ,z=0,w)  -  RSS  A,  +  RDS  A2  +  RDD  A3  (A5.1) 

d4(r,zM ),w)  -  TSS  A4  +  TDS  Ah  , 

where  dt  is  the  vertical  displacement  (positive  upward),  dr 
is  the  radial  displacement,  and  is  the  tangential  dis¬ 
placement  (positive  in  a  direction  clockwise  from  north). 
The  functions  ZSS,  ZDS,  ZDD,  RSS,  RDS,  RDD,  TSS, 
and  TDS  together  with  ZEP  and  REP  are  the  ten  Green’s 
functions  required  to  calculate  a  wave  field  due  to  an  arbi¬ 
trary  point  dislocation  source  or  point  explosion  buried  in 
a  plane  layered  medium  (Wang  and  Herrmann,  1980; 
Herrmann  and  Wang,  1985).  As  before,  let  u=(uJttif ,uz) 
and  be  the  dislocation  vector  and  vector 

normal  to  the  fault  plane,  respectively.  Note  that  (15)  and 
(16)  are  identical  to  the  formulation  used  by  Herrmann 
and  Wang  (1985),  where  our  u  equals  their  f  and  our  v 
equals  their  n  (1  =  x-axis,  2  =  y-axis,  3  =  z-axis).  Then 

A  j— (u,ut  —  cos(2az)+(u1i^+u>i'I)  sin(2az) 

A2=(uIi',+u,i'1)  cos(az)+(uy  L'J+uJvt )  sin(az) 

A3—utut  (A5.2) 

A4=(tiIt/I—u>i/y)sin(2az)—(uJl/t+v)i'J)  cos(2az) 

"i+“,l'x)  sin(«z)-(u,i/,-Hi1i/,)  cos(az)  , 

where  az  is  the  azimuth  of  observation.  Equivalently, 

A  ,=  y(A/IZ  — A/,v  )  cos(2az)+A/zll  sin(2az) 

A  2»=  A/IZ  cos(az)+A/JZ  sin  (az) 

A  3  (A5.3) 

A  A—  ~{MIZ  — Afj, )  sin(2az  )— A/ZJ  cos(2 az) 

A5=  —  A/>z  cos(az)+Aflz  sin(az)  . 

These  equations  are  identical  to  (A5.2)  which  can  be 
proven  by  using  (18)  together  with  (15)  and  (16).  Note 
that  the  coefficients  given  in  (A5.3)  agree  with  the 
moment  tensor  elements  as  defined  by  Aki  and  Richards 
(1980;  (A5.3)  differs  in  sign  with  the  coefficients  of  Langs¬ 
ton  (1981)  due  to  conventions  on  displacements  and 
Green’s  functions). 

Note  that  either  definition  of  the  coefficients  of  the 
Green’s  functions  can  be  used  for  the  calculation  of  the 
displacement  at  the  free  surface,  depending  on  whether 
the  focal  mechanism  or  the  moment  tensor  is  given.  Here, 
equations  (A5.3)  and  (A5.1)  are  used  in  order  to  develop 
an  inversion  scheme  for  the  moment  tensor  elements.  We 


regroup  and  assume  the  presence  of  an  isotropic  com¬ 
ponent  (ZEP  #  0,  REP  ^  0): 


rfz(r,z=0,u/)  -  A/„ 


+  M, 


» 


+  A/, 

+ 

+  M, 

+ 


ZSS  ,  ZDD  .  ZEP 

—  cos(2az) - _  + 

-ZSS  ,  ZDD  ,  ZEP 

__cos(2(u) 

ZEP 


|zSS  sin(2az)j 
| ZDS  cos(az)J 
Mft  |zZ)S  sin(az)j 


(A5.4) 


dr(r,z«0,w)  -  Ma 

RSS  , 

—z—  cos(2oz)  - 

+  Mn 

—RSS 

cos(2az) 

2 

+  Mu 

REP 

3 

+  M„ 

RSS  sin(2az)  j 

+ 

RDS  cos( az)  j 

+  Mt> 

RDS  sin(az)J 

d+(r  ,z=0,w)  =  A/„ 

TSS  .  , 

— - —  sin(2az) 

+  Mjj 

-TSS 

sin(2az) 

2 

+  Mt) 

to 

£ 

1 

cos(2az )  j 

+  M„ 

TDS  sin(az)j 

+  Mt, 

-TDS 

cos(az)  J 

RDD  REP 
2  3  J 

RDD  REP 
2  3 


(A5.5) 


(A5.6) 


Equations  (A5.4),  (A5.5),  and  (A5.6)  each  set  up  a 
moment  tensor  inversion  scheme.  Equations  (A5.4)  and 
(A5.5)  are  formulated  for  the  general  case  where  the  inver¬ 
sion  expects  a  moment  tensor  that  is  a  composition  of  an 
isotropic  part  and  a  deviatoric  part.  An  inversion  based 
on  transverse  data,  (A5.6),  cannot  resolve  A/„.  In  such  a 
case,  we  assume  that  the  moment  tensor  is  purely  devia¬ 
toric  and  constrain  Mtl  «=  -  (Mu  -I-  Mn).  The  same  con¬ 
straint  can  be  applied  to  (A5.4)  and  (A5.5)  in  the  case  of 
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an  inversion  that  looks  for  a  pure  deviatoric  moment  ten¬ 
sor  (set  formally  ZEP  =  REP  =  0  in  (ASA)  and  (A5.5), 
Dziewonski  el  al  . ,  1081). 

From  the  last  three  equations  we  see  that  the 
observed  displacement  at  the  free  surface  is  a  linear  com¬ 
bination  of  the  station  specific  Green's  functions,  within 
the  square  brackets,  with  the  moment  tensor  elements  as 
scalar  multipliers.  We  also  note  that  if  the  source  time 
function  is  known  and  a  point  source  approximation  is 
acceptable,  the  moment  tensor  elements  are  independent 
of  frequency  (linear  inversion)  and  similar  equations  arise 
relating  observed  time  histories  to  temporal  Green’s  func¬ 
tions  within  the  square  brackets. 

Next,  we  performed  a  simple  moment  tensor  inver¬ 
sion  using  the  vertical  component  of  synthetic  teleseismic 
P-wave  first  motion  peak  amplitudes  as  suggested  by 
Stump  and  Johnson  (1977).  We  assumed  a  pure  devia¬ 
toric  source  (ZEP  =  0  in  (A5.4)). 

Let  az,,  ...,  az„  be  azimuths  of  n  different  stations. 
Then  the  expressions  in  the  square  brackets  of  (A5.4) 
define  components  of  a  matrix  as  a^az,)  ,  ...,  dlS(az,)  for 
the  i-th  azimuth.  A  system  of  linear  equation  arises: 

<hs(«2i)  F,w  ' 

JJZZ 

Afyy 

-  A/Iy  (A5.7) 

K, 

n)  “„l(«Al)  “n5(“OJ  [A/yc  . 

For  observations  at  more  than  5  distinct  azimuths, 
the  system  (A5.7)  is  overdetermined.  The  solution  can  be 
reached  by  the  classical  least  squares  approach.  The  five 
moment  tensor  elements  can  be  determined  by  using  the 
numerical  stable  singular  value  decomposition.  We 
imposed  the  deviatoric  constraint  Mu  =  -  (Af„  +  A/vy ). 
Hence  the  inversion  gives  a  purely  deviatoric  source. 
However,  we  were  not  constraining  one  eigenvalue  as  zero 
(double  couple),  letting  the  inversion  tell  us  about  double 
couple  and  CLVD  components.  The  eigenvalues  and 
eigenvectors  can  be  calculated  using  the  Householder 
transformation  with  further  QL  decomposition.  The 
implementation  of  these  numerical  concepts  was  done 
using  code  by  Press  el  al  (1987). 

In  the  following,  some  results  of  inverting  synthetic 
data  are  presented.  First,  Green's  functions  were  calcu¬ 
lated  using  a  Haskell  formalism  for  a  simple  half-space 
model  (  Vp  =  8  km/sec  ,  V,  =  4.6  km/sec  and  p  =  3.3 
g /cm3,  h  =  30  km  ).  Figure  A. 3  shows  the  three  basic 
Green’s  functions  ZSS,  ZDD,  and  ZDS.  The  assumed 
focal  mechanism  (Figure  A. 4:  strike  =  180”,  dip  =  40°, 
slip  =  110°)  is  the  same  as  used  by  Herrmann  (1975,  Fig¬ 
ure  2).  Teleseismic  P-wave  first  motions  were  synthesized 
at  12  equidistant  azimuths  (epicentra!  distance  =  50°). 
Note  that  an  instrument  response  was  not  included  in  the 
synthetics.  Due  to  the  simple  model  and  the  fact  that  all 
stations  are  equidistant  from  the  source,  a  correction  for 
anelastic  attenuation  (  ( *  =  0.7  )  or  geometrical  spread¬ 
ing  is  not  required.  A  correction  for  an  extended  source  is 
not  necessary  since  the  moment  used  is  1020  dyne-em  and 
the  duration  of  the  source  time  function  is  0.2  sec.  We 
used  (A5.4)  for  time  domain  measurements. 


Fig.  A. 3.  Synthetic  Green’s  functions  ZSS,  ZDD.  and  ZDS 
(Herrmann  and  Wang,  1985)  for  a  half-space  model 
(  Vf  —  8  km/sec,  V,  =  4.6  km/sec,  p  =  3.3 
g /cm3,  h  =  30  xm,  t'  =  0.7)  calculated  by  using 
the  Haskcli  formalism.  The  time  window  ranges 
from  4.0  to  55.1  sec  (dt  =  0.05  sec).  Maximum 
amplitudes  are  in  cm  (Appendix  V). 


N 


Fig.  A.1.  Assumed  focal  mechanism  for  the  synthetic 
seismograms:  strike  =  180°,  dip  =  40°,  slip  = 
110°  (Appendix  V). 
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TABLE  A.7:  RESULTS  OF  THE  MOMENT 
TENSOR  INVERSION  (MAJOR  COUPLE) 


Case  0 

Case  ! 

Case  I! 

Case  III 

Case  IV 

CaseV 

0 

-0.037 

-0  050 

-0.109 

-0.202 

0.301 

M„ 

-0.925 

-0  902 

-0  951 

-0.966 

-1.023 

-1.091 

A/„ 

0.925 

0  939 

1.002 

1  075 

1  225 

0  791 

M„ 

-0.220 

-0.199 

-0.200 

-0.194 

-0  176 

0.257 

M„ 

-0.262 

-0  262 

-0.260 

-0.264 

-0.257 

-0  172 

M„ 

-0.163 

-0  168 

-0.162 

-0.168 

-0.156 

-0  324 

EV(NULL) 

0.00 

-0.04 

-0  05 

-Oil 

-0.20 

0.26 

EV(T) 

1  00 

1.01 

1.07 

1  14 

1.28 

0.92 

EV(P) 

-1.00 

-0.97 

-1.02 

-1.03 

-1.08 

-1.18 

%of  DC 

100 

92 

90 

81 

69 

56 

%  of  CLVD 

0 

8 

10 

19 

31 

44 

A/0 

1.00 

0  99 

1.04 

1  09 

1  19 

1.07 

STRIKE 

180.0 

179.5 

179  2 

177  8 

176  4 

211  8 

DIP 

40.0 

39  7 

40.1 

39  9 

40  5 

40  8 

SLIP 

110  0 

109  0 

107  8 

106.1 

103  2 

123.1 

STRIKE 

331  6 

335  4 

336.5 

337.2 

339  3 

351.1 

DIP 

52.8 

52.9 

52.2 

51.9 

508 

56  8 

SLIP 

74.0 

74.9 

75.8 

77.0 

70.0 

64  8 

T  (TRD) 

192.7 

194.9 

194.7 

196.8 

198  2 

209  7 

T  (PLC) 

75.8 

76.2 

77.1 

78.1 

800 

67.3 

P  (TRD) 

75  9 

76.1 

76.7 

76.4 

77  1 

98  8 

P  (PLG) 

6.6 

6.7 

6  2 

6  1 

5.2 

8  5 

For  Case  0,  moment  tensor  elements  are  calculated  from  (18)  assuming  a  double  couple  source 
(strike  =  180',  dip  —  40'.  slip  =  1 10').  The  eigenvalues  of  the  moment  tensor  corresponding  to 
the  null-,  T-,  and  P-axes  are  shown  as  EV(NULL),  EV(T).  and  EV(P),  respectively.  Equation  (38) 
is  used  to  determine  the  percentage  of  double  couple  or  CLVD  from  the  eigenvalues  of  the  moment 
tensor  The  seismic  moment  is  calculated  using  (20)  The  orientation  of  the  fault  plane  and  auxili¬ 
ary  plane  is  given  together  with  the  trend  and  plunge  of  the  T-  and  P-axes  (Herrmann,  1975) 
Cases  ]  -  IV  are  for  additive  pseudo-random  noise  (  0  %,  14  %,  28  %,  and  56  %,  respectively)  in 
the  synthetic  seismograms  at  12  different  azimuths  Case  V  assumes  that  one  of  the  12  seismo¬ 
grams  has  a  reversed  polarity  (0  %  pseudo  random  noise). 

Table  A.7  displays  the  inversion  results  for  the  major 
couple.  The  moment  tensor  elements,  the  percentage  of 
double  couple  and  CLVD,  the  seismic  moment,  and  the 
focal  mechanism  parameters  are  shown.  For  Case  0,  the 
moment  tensor  elements  were  calculated  from  the  given 
fault  plane  solution  and  (18).  Next,  three  experiments 
were  performed:  1.)  synthetic  seismograms  were  calculated 
using  the  Haskell  method  (Case  I).  Figure  A.5a  shows  the 
vertical  component  of  a  synthetic  seismogram  at  an 
azimuth  of  0  degrees.  2.)  Different  amounts  of  pseudo¬ 
random  noise  were  added  to  the  synthetic  seismograms 
calculated  in  Case  1  with  amplitudes  of  ±  0.25  X  10-9  cm 
(Case  II,  Figure  A. 5b),  ±  0.5  X  10-9  cm  (Case  III,  Figure 
A. 5c),  and  ±  1.0  X  10-9  cm  (Case  IV,  Figure  A.5d). 
Averaged  over  the  12  azimuths,  these  noise  levels 
correspond  to  14  %,  28  %,  and  56  %  pseudo-random 
additive  noise,  respectively.  3.)  The  final  experiment 
(Case  V)  relates  to  possible  polarity  errors  of  seismo¬ 


graphs.  Hence  it  was  assumed  that  one  of  the  12  seismo¬ 
grams  of  Case  I  had  a  wrong  polarity. 

The  theoretical  focal  parameters  (Case  0)  agree 
within  the  measurement  errors  with  the  observed  ones 
(Case  I).  This  justifies  the  technique.  The  effect  of  noise  is 
to  severely  distort  the  moment  tensor  elements.  The  iso¬ 
tropic  moment  tensor  components  seem  to  be  more  sensi¬ 
tive  to  noise  than  the  deviatoric  ones.  Notice  that  the 
moment  tensor  gains  a  contribution  of  a  CLVD  due  to  the 
noise.  The  percentage  of  CLVD  versus  double  couple 
increases  with  increasing  noise.  The  effect  of  random 
noise  on  the  fault  plane  solution  that  is  derived  from  the 
moment  tensor  elements  is  minor;  i.e.  the  fault  plane 
solution  for  the  major  double  couple  is  very  close  to  the 
original  focal  mechanism.  However,  with  increasing  noise, 
the  fault  plane  solution  deteriorates.  8  %  polarization 
errors  in  otherwise  perfect  data  lead  to  worse  results  than 
56  %  additive  random  noise  (Case  IV)-  A  doubling  of  the 
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2. 2&7E-09 


2.981tE-09  d 


Fig.  A.5.  Vertical  components  of  synthetic  teleseismic 
seismograms  at  50  degrees  and  azimuth  of  0°.  The 
time  window  ranges  from  4.0  to  20.6  sec  (dt  = 
0.05  sec).  Maximum  amplitudes  are  in  cm.  (a)  No 
pseudo-random  noise  added;  (b)  -  (d)  pseudo¬ 
random  noise  is  added  with  amplitudes  of 
±0.25  X  10'®  cm,  ±0.50  X  10"®  cm,  ±1.0  X  10“® 
cm,  respectively  (Appendix  V). 


polarization  errors  gives  meaningless  results.  Due  to  the 
setup  of  the  experiment,  a  minor  couple  would  be  a  pure 
artifact  of  the  noise. 
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Ground  roll:  rejection  using  adaptive  phase  matched  filters 
Robert  B.  Herrmann*  and  David  R.  Russell** 

ABSTRACT 

The  technique  of  phase  match  filtering  dispersive  surface  waves  is  extended  to 
permit  an  adaptive,  iterative  process  by  which  the  signal  itself  in  a  seismic  trace 
designs  a  filter  to  remove  the  surface  wave.  The  technique  is  robust  and  well  behaved, 
and  requires  the  specification  of  only  simple  parameters  for  its  operation. 

The  technique  is  applied  to  data  sets  from  three  regions,  representing  a  wide 
range  in  the  ratio  of  surface-wave  noise  to  exploration  signal.  The  technique  works 
very  well  with  poor  data  sets  and  also  improves  good  data  sets.  Since  the  technique  is 
applied  to  individual  traces,  it  works  in  situations  for  which  f-k  filtering  might  not  be 
feasible  due  to  poor  spatial  sampling.  The  technique  is  computationally  more  intensive 
than  recursive  digital  bandpass  filtering  of  individual  traces,  but  is  less  intensive  than 
filtering  in  the  f-k  domain. 
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INTRODUCTION 

Surface-wave  noise,  the  ground  roll,  on  reflection  seismograms  is  always  a  part  of 
land  based  data  acquisition.  As  such  it  has  always  been  something  to  be  eliminated  in 
order  to  focus  on  the  underlying  reflection  data.  Dobrin  (1951)  studied  the  ground  roll, 
showing  that  it  possessed  the  properties  of  dispersive  surface  waves.  Al-Husseini  et  al 
(1981)  inverted  the  surface-wave  data  obtained  in  special  tests  to  characterize  the 
ground  roll  in  a  region  of  eastern  Saudi  Arabia  and  thus  to  provide  the  necessary  infor¬ 
mation  for  proper  group  array  design  for  its  elimination.  Under  optimal  conditions, 
proper  group  array  design  succeeds  in  attenuating  the  low  phase  velocity  surface-wave 
arrivals,  enhancing  the  high  phase  velocity  reflections. 

Earthquake  seismology  has  historically  taken  the  opposite  approach,  that  the  sur¬ 
face  wave  is  a  very  useful  signal  for  defining  earth  structure  as  well  as  the  seismic 
source.  Because  of  this  emphasis  many  techniques  for  surface-wave  analysis  have 
been  developed  in  this  field.  Our  technique  makes  use  of  phase  matched  filtering  of 
surface  waves  (Herrin  and  Goforth,  1977;  Goforth  and  Herrin,  1979;  and  Beresford- 
Smith  and  Mason,  1980).  In  earthquake  seismology,  the  objective  of  phase  matched 
filtering  is  to  isolate  a  dispersed  surface  wave  mode  from  a  noise  background  consist¬ 
ing  of  body-wave  arrivals  or  secondary  surface-wave  arrivals  due  to  multipathing.  The 
isolated  surface  wave  has  a  smoother  amplitude  spectrum  and  also  a  better  defined 
dispersion  for  use  in  source  and  earth  structure  studies,  respectively.  From  this  point  of 
view,  the  typical  exploration  reflection  signal  is  viewed  as  noise.  The  essence  of  the 
technique  described  in  this  paper  is  the  recognition  that  the  surface-wave  signal  in 
earthquake  seismology  is  the  exploration  noise,  and  vice  versa. 

Phase  matched  filtering 

Consider  a  frequency  domain  representation,  S(f ),  of  a  time  domain  signal,  s(t) 
at  a  distance  x ,  consisting  of  N  arrivals: 

S(f)='ZAi(f)e~jki(S)x.  (1) 

i=l 

Here,  a  propagating  wave  notation  is  used  to  represent  each  arrival.  At  (f )  and  kt  (f ) 
are  the  complex  amplitude  and  wavenumber  spectra  of  the  i’th  arrival.  If  the 
wavenumber  kt(f)  is  purely  linear  in  frequency,  then  its  corresponding  time  domain 
signal  is  just  time  shifted  as  a  function  of  distance;  otherwise  it  is  dispersed.  In  addi¬ 
tion,  as  long  as  we  focus  on  a  single  distance,  equation  (1)  is  general  enough  to  also 
include  a  non-propagating  noise  components. 

Suppose  further  that  the  m'th  arrival  is  not  desired,  and  that  the  wavenumber 
function  Km(/)  is  a  reasonable  approximation  to  the  unknown  km(f).  Now  multiply 
S  (f )  by  an  inverse  propagating  wave  function  to  yield 

S(f)ejKm(f)x  =  £  Ai(f)e~j[ki(fhKm(/)]x  (2) 

+  Am(f)e-j[k”(fhK~(f)]x. 

If  the  terms  [k,-  (f  )-Km  (f  )]x  for  i*m  are  sufficiently  different  from  zero,  and  if  the 
term  [km(f  )-Km(f  )]x  is  sufficiently  close  to  zero,  then  the  inverse  transform  of 
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equation  (2)  will  yield  a  pulse  at  zero  lag,  which  will  be  the  compressed,  or 
undispersed,  m'th  arrival.  The  other  arrivals  will  be  either  spread  out,  or  if 
compressed,  not  compressed  as  much,  and  certainly  not  located  at  zero  lag. 

In  earthquake  seismology,  the  signal  so  compressed  and  shifted  to  zero  lag  is  the 
one  of  interest.  It  is  windowed  about  zero  lag,  transformed  into  the  frequency  domain, 
uncompressed  using  the  dispersion  operator 

(3) 

and  then  inverse  transformed  into  the  time  domain.  This  operation  yields  a  clean 
surface-wave  signal  (Herrin  and  Goforth,  1977;  Goforth  and  Herrin,  1979).  To  obtain 
the  exploration  signal  one  can  either  subtract  the  isolated  m'th  arrival  from  the  original 
time  series,  the  technique  used  here,  or  mute  the  inverse  transform  of  equation  (2) 
about  zero  lag,  and  then  take  a  Fourier  transform  of  the  resultant  trace,  apply  the 
operator  of  equation  (3),  and  then  inverse  transform  the  result  to  the  time  domain 
(Beresford-Smith  and  Rango,  1988;  Saatpilar  and  Canitez,  1988).  The  result  in  either 
case  will  now  be  a  multi-arrival  signal  lacking  the  m'th  arrival. 

While  conceptually  very  simple,  the  success  of  this  technique  lies  in  the  correct 
specification  of  the  function  Km  (f ).  If  multitrace  data  are  available,  and  if  a  coherent 
k-f  or  /7-co  (McMechan  and  Yedlin,  1981)  display  can  be  made,  then  this  may  be 
possible.  The  resulting  k m(f)  would  represent  a  spatial  average  over  the  data  set, 
which  may  not  be  appropriate  if  there  are  lateral  variations  in  the  surface-wave  disper¬ 
sion.  On  the  other  hand,  if  the  objective  is  only  to  remove  the  undesired  surface-wave 
signal  and  not  to  define  the  correct  phase  velocity  dispersion  function,  then  this  same 
technique  can  be  applied  using  single  trace  data,  e.g.,  treating  each  trace  as  an 
independent  data  set.  The  trick  is  to  make  the  processing  iterative,  allowing  the  data 
itself  to  define  its  own  dispersion  operator  by  improving  an  initial  estimate.  This  is  in 
fact  the  heart  of  the  Herrin  and  Goforth  (1977)  and  Goforth  and  Herrin  (1979)  tech¬ 
nique. 

In  order  to  reject  a  signal,  it  must  be  isolated.  Compressing  it  perfectly  should 
yield  a  zero  phase  wavelet  centered  at  zero  lag.  The  degree  to  which  this  is  not  true, 
indicates  the  subtle  differences  between  the  unknown  km(f)  and  the  trial  function 
k m(f).  By  windowing  the  compressed  signal  about  zero  lag,  Fourier  transforming  it, 
and  unwrapping  the  phase  delay,  the  difference  can  be  defined.  This  phase  difference 
is  then  used  to  adjust  the  k m(f). 

To  control  this  procedure,  we  force  the  function  k m(f)  to  be  initially  a  simple 
dispersion  function  by  approximating  it  by  B-splines  with  a  small  number  of  nodes.  At 
each  iteration,  a  least-squares  B-spline  is  fit  to  the  phase  differences,  to  yield  a  simple, 
smooth  perturbation  on  the  current  estimate  of  Km  (f ).  This  smoothing  technique  main¬ 
tains  the  assumption  that  a  smooth  dispersion  operator  is  required  to  model  surface 
waves. 

To  start  the  process,  an  initial  estimate  of  Km(f)  is  required.  For  simplicity  of 
use,  we  use  the  following  approach.  By  visual  examination  of  all  traces  or  by  applica¬ 
tion  of  multiple  filter  techniques  (Dziewonski  et  al.,  1969;  Herrmann,  1973),  an  esti¬ 
mate  is  made  of  the  group  velocity  window  and  bandwidth  associated  with  each  signal 
that  is  to  be  removed.  Let  the  group  velocity  window  be  defined  by  the  group 
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velocities  Uq  and  U  j,  and  let  the  corresponding  frequencies  associated  with  these  velo¬ 
cities  be/o  and  /  lt  respectively.  Following  Saatpilar  and  Canitez  (1988),  we  approxi¬ 
mate  the  group  velocity  dispersion  function  by  requiring  that  the  group  slowness  be  a 
linear  function  of  frequency.  Since  the  relation  of  wavenumber  k  (f )  to  group  slowness 


CJ  1(/)  is  given  by  the  definition  k(f)  =  k0  +  2k  ju  1  Qt )dx ,  the  initial  k m(f)  esti- 

/o 

mate  is  defined  once  we  specify  a  k0.  We  choose  the  value 


kn  =  2k 


center  o: 


equation  (2). 


fo+f] 


UK1  +  ur1 


to  guarantee  that  the  desired  signal  at  the  effective 
the  group  velocity  window  is  shifted  to  zero  lag  on  the  first  application  of 


In  applying  this  technique,  the  input  trace  is  windowed  within  the  specified  group 
velocity  window  by  using  a  window  with  a  width  twice  the  longest  period  in  order  to 
focus  the  processing  on  the  desired  signal,  and  three  iterations  are  performed.  The 
number  of  Fast  Fourier  Transforms  required  are  2(1  +  Nj),  where  Nj  is  the  number  of 
iterations. 


This  processing  technique  is  adaptive  since  the  surface  wave  signal  itself 
improves  the  dispersion  function.  The  technique  is  robust  since  an  initial  poor  guess  of 
the  dispersion  curve  may  yield  nothing  of  consequence  near  zero  lag  when  the  signal 
is  reduced  to  zero  distance,  in  which  case  the  phase  match  filtered  signal  is  almost 
non-existent,  and  the  result  of  the  signal  subtraction  is  essentially  the  same  as  the  ini¬ 
tial  signal. 


DATA  PROCESSING 

Oklahoma 

As  part  of  a  sensor  evaluation  test,  four  data  sets  were  collected  at  a  field  test  site 
at  Oklahoma.  Both  dynamite  and  vertical  vibrator  sources  were  used  with  receiver 
group  arrays  spaced  15.2  m  apart  between  152.4  and  441.9  m  from  the  source.  Each 
group  array  consisted  of  six  geophones.  The  vibrator  acted  at  the  surface  while  the 
dynamite  was  buried  at  a  depth  of  15.2  m.  The  four  vertical  component  data  sets 
acquired  are  as  follow: 


DynGrp  -  dynamite  source  with  a  30.5  m  geophone  group  array 
DynNoGrp  -  dynamite  source  with  a  1  m  geophone  group  array 
VibGrp  -  vibrator  source  with  a  30.5  m  geophone  group  array 
VibNoGrp  -  vibrator  source  with  aim  geophone  group  array 

The  receiver  group  arrays  were  not  designed  to  maximally  reduce  the  surface  waves  ; 
at  most  they  reduce  the  surface-wave  signal  by  6-10  db  over  the  usable  signal  range 
(5-85  Hz). 

Figure  1  presents  the  input  time  histories,  showing  1000  ms  of  data  sampled  at  4 
ms  intervals.  Within  each  panel,  the  trace  on  the  left  is  at  a  distance  of  152.4  m  and 
the  one  at  the  right  is  at  a  distance  of  441.9  m.  An  automatic  gain  correction  with  a 
500  ms  window  has  been  applied  to  each  plotted  trace.  In  this  figure,  a  set  of 
reflections  is  seen  at  a  two  way  traveltime  of  about  500  ms.  This  is  from  the  well 
known  Woodford  shale  formation  at  about  610  m.  In  addition  well  developed  surface 
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waves  are  observed.  The  data  sets  are  arranged  such  that  the  ratio  of  reflection  signal 
to  surface  wave  noise  is  greatest  on  the  far  left  panel  and  decreases  to  the  right.  This 
is  well  understood  theoretically  in  that  buried  sources  excite  the  surface  wave  less  than 
surface  sources,  and  receiver  group  arrays  of  reasonable  lateral  dimensions  also  reduce 
the  low  phase  velocity  arrivals.  We  note  that  the  500  ms  reflection  is  not  apparent  in 
the  VibNoGrp,  Figure  Id,  data  set  at  all. 

One  technique  for  reducing  the  surface  wave  in  a  record  uses  the  fact  that  the 
surface  wave  usually  has  a  lower  frequency  content  than  the  reflections.  In  this  area 
the  peak  surface-wave  amplitude  is  at  IS  Hz.  Application  of  a  high  pass  Alter  to  the 
record  section  will  attenuate  the  surface  wave.  Figure  2  shows  the  result  of  applying  a 
2  pole,  zero-phase  high-pass  Butterworth  Alter  with  a  low-cut  comer  frequency  at  32 
Hz  to  the  Aeld  data.  It  is  obvious  that  the  low  frequency  surface  wave  is  virtually 
eliminated  in  the  case  of  the  data  set  with  the  best  S/N  ratio,  DynGrp,  Figure  2a. 
However,  low  velocity  surface-wave  and  shear-wave  arrivals  can  be  seen  in  the  other 
time  histories.  The  comer  frequency  of  the  Alter  could  be  increased  to  further  reduce 
the  surface-wave  arrivals,  but  this  would  in  turn  reduce  the  bandwidth  of  the 
reflections. 

To  test  the  effectiveness  of  the  phase  matched  rejection  Alter  technique,  we  pro¬ 
cessed  each  trace  through  three  different  Alters,  with  control  parameters  given  in  Table 
1.  For  each  filter  so  defined,  three  iterations  were  used  to  adaptively  refine  the  filter, 
phase  delays  were  smoothed  using  a  5  point  least  square  B-spline.  An  AGC  was 
applied  to  the  data  prior  to  phase  match  filtering  each  trace  to  amplify  the  signal  at  the 
end  of  the  trace.  The  AGC  gain  as  a  function  of  time  was  saved  so  that,  after  the 
phase  match  filtering,  the  effects  of  the  AGC  could  be  removed  to  preserve  the  true 
amplitudes  of  the  underlying  reflections. 

Figures  3  and  4  show  the  initial  data  sets  and  the  output  of  each  phase  of  pro¬ 
cessing  for  the  best,  DynGrp,  and  worst,  VibNoGrp,  data  sets  in  terns  of  S/N, 
respectively.  The  parameters  of  the  first  pass  were  determined  following  the  applica¬ 
tion  of  standard  techniques  for  determining  group  velocities  from  single  traces  ( 
Dziewonski  et  al,  1969)  The  first  pass  clearly  removes  the  fundamental  mode  surface 
wave.  For  the  DynGrp  data  set,  further  processing  does  not  change  the  character  of 
the  record  section  in  Figure  3. 

The  processing  parameters  used  in  the  second  pass  were  an  attempt  to  remove  the 
direct  shear-wave  arrivals  observed  between  the  first  P  break  and  the  500  ms  reflection 
that  is  prominent  in  the  DynNoGrp,  VIbGrp  and  VibNoGrp  data  sets.  This  has  a 
higher  frequency  content  than  the  fundamental  mode  surface  wave  and  travels  at  a 
higher  phase  velocity.  The  result  of  this  pass  is  quite  significant  in  the  case  of  Vib¬ 
NoGrp  for  which  the  500  ms  reflection  now  emerges.  The  last  pass  is  an  attempt  to 
remove  the  very  low  group  velocity  arrivals  that  may  be  generated  by  the  air  wave. 

Figure  5  compares  the  results  of  applying  the  three  pass  phase  matched  filtering 
operation  to  the  original  data  sets.  Comparing  this  figure  to  initial  data  sets  shown  in 
Figure  1,  it  is  obvious  that  there  has  been  a  significant  improvement  in  the  reflection 
signal  to  surface-wave  noise  ratio. 

To  see  how  these  results  compare  to  that  of  simple  high-pass  filtering,  Figure  6 
shows  the  result  of  high-pass  filtering  the  phase  matched  processing  output  using  the 
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same  filter  used  in  Figure  2.  Comparing  this  to  Figure  2  shows  the  value  of  reducing 
the  shear  wave  arrivals  through  phase  match  filtering.  In  addition  the  faster  surface 
wave  seen  in  the  DynNoGrp,  VibGrp  and  VibNoGrp  data  sets,  Figures  2b,  2c,  and 
2d,  respectively,  has  been  eliminated. 

This  data  set  demonstrates  how  phase  matched  filtering  can  be  applied  to  high 
resolution  data  collection  strategies  using  either  small  or  zero  length  group  arrays.  A 
more  positive  view  is  that  good  results  may  be  obtained  using  fewer  instruments  in  the 
group  arrays  as  long  as  there  is  numerically  resolvable  signal  beneath  the  surface 
wave. 

West  Texas 

The  second  data  set  comes  from  a  high  density  mini-refraction  profile  in  west 
Texas  which  was  notable  for  its  lack  of  reflections,  but  which  was  useful  for  other  stu¬ 
dies  because  of  the  sharp  P-wave  first  breaks.  The  displays  of  the  data  set  and  phase 
matched  outputs  are  AGC’d  with  a  500  ms  window.  Figure  7a  presents  the  original 
data  set.  Traces  are  at  15.2  m  intervals  from  15.2  to  731.5  m.  Each  trace  consists  of 
1000  samples  at  a  2  ms  sampling  interval.  The  initial  P  waves  and  their  reverberations 
near  the  surface  are  apparent,  as  is  a  high  frequency  air  wave  arrival,  e.g.,  at  2.00  s  at 
a  distance  of  701.0  m.  In  addition  there  are  a  number  of  overlapping  dispersed  low  fre¬ 
quency  arrivals. 

The  phase  matched  filter  was  run  as  for  the  Oklahoma  data  set  except  the  disper¬ 
sion  parameters  of  Table  2  were  used.  The  parameters  for  the  first  pass  were  chosen 
on  the  basis  of  a  /?-/  stack  and  group  velocity  analysis.  The  next  passes  were 
designed  to  remove  the  surface  wave  following  the  air  wave  and  the  direct  shear  wave. 
The  final  pass  was  an  attempt  to  remove  the  dispersed  reverberation  following  the 
direct  P  arrival.  The  reduction  in  surface-wave  signal  is  evidenced  by  the  relative 
enhancement  of  the  air  wave  arrival  in  the  AGC  display.  The  interesting  aspect  of  this 
data  set  was  the  dispersed  surface  wave  following  the  air  wave.  These  arrivals  were 
set  up  by  a  moving  surface  source,  the  air  wave  propagating  across  the  array  at  300 
m/s  (Press  and  Ewing,  1951;  Mooney  and  Kaasa,  1962;  Ren6  et  al,  1986).  The  third 
pass  chose  a  group  velocity  window  from  the  air  wave  to  a  point  later  in  the  trace. 
The  final  result,  Figure  7b  demonstrates  a  significant  reduction  in  the  surface  wave, 
whether  set  up  directly  by  the  source  or  indirectly  by  the  propagation  air  wave  distur¬ 
bance. 

Permafrost 

Data  sets  in  the  arctic  regions  are  notorious  for  poor  signal  to  noise,  and  often  are 
excellent  candidates  for  phase  match  filters.  Barton  et  al  (1986),  McConnell  et  al 
(1986),  and  Beresford-Smith  and  Rango  (1988)  discuss  the  problems  with  data 
acquisition  on  floating  ice  sheets.  Our  data  set,  Figure  8a,  was  acquired  on  land,  and 
is  not  overwhelmed  by  the  dispersive  flexural  waves.  The  traces  range  in  distance  from 
1542  m  on  the  left  to  100  m  on  the  right.  There  are  some  very  strong  undispersed 
arrivals  with  group  velocities  of  1400  m/s  and  300  m/s.  The  large  amplitude  of  these 
arrivals  is  indicated  by  the  apparent  muting  introduced  immediately  prior  because  of 
the  AGC  process.  These  signals  are  also  very  narrow  band  with  center  frequency  near 
15  Hz.  The  result  of  processing  this  data  set  with  the  parameters  of  Table  3  is  shown 
in  Figure  8b.  The  air  wave  is  still  present  as  a  very  high  frequency  arrival,  but  the 
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1400  m/s  arrival  is  virtually  eliminated.  The  reduction  of  peak  amplitude  is  indicated 
by  the  lack  of  a  muting  by  the  AGC. 

This  data  set  demonstrates  the  inherent  control  of  our  technique,  in  that  it  can  be 
applied  to  a  very  narrow  group  velocity  window  and  that  if  no  dispersion  is  present,  it 
defaults  to  bandlimited  mute. 


CONCLUSIONS 

A  technique  well  known  in  earthquake  seismology  has  been  applied  to  an 
exploration  problem  with  the  underlying  idea  that  the  earthquake  noise  is  the  explora¬ 
tion  signal  and  that  the  earthquake  signal  is  the  exploration  noise.  The  adaptive  phase 
matched  filter  technique  differs  from  k-f  filter  techniques  in  that  each  trace  is 
independently  analyzed.  This  means  that  our  method  will  work  for  data  sets  which 
would  be  severely  aliased  in  k  -f  analysis.  To  some  extent  it  is  similar  to  time  variant 
spectral  whitening  technique  described  by  Yilmaz  (1987)  in  that  undesired  signal  is 
removed  from  a  single  trace  in  a  time  variant  manner.  However,  our  technique  essen¬ 
tially  subtracts  a  coherent  noise  from  the  trace,  leaving  the  underlying  signal  undis¬ 
torted  in  frequency  and  time,  thus  preserving  true  amplitude. 

A  similar  technique  has  been  applied  to  the  problem  of  seismic  data  acquisition 
on  ice  to  remove  the  dispersive  flexural  wave  (Barton  et  al,  1986;  McConnell  et  al, 
1986;  Beresford-Smith  and  Rango,  1988).  These  approaches  used  a  phase  matched 
filter  and  used  an  f-k  analysis  to  define  the  dispersion  operator,  though  Barton  et  al 
(1986)  recognized  the  problems  with  lateral  changes  in  dispersion.  They  then  applied 
a  phase  matched  rejection  filter  to  the  data  sets.  The  adaptive  phase  matched  filter 
technique  presented  here  does  not  require  such  a  well  defined  dispersion  operator. 
Only  a  reasonable  estimate  is  required.  We  used  simple  group  velocity  windows, 
which  are  easily  chosen  by  an  analyst.  This  simple  approach  works  since  the  dispersed 
seismic  signal  itself  defines  its  own  rejection  filter. 

The  processing  shown  was  done  on  a  number  of  different  UNIX  (TM  AT&T 
Technologies)  machines.  Using  the  present  unoptimized  code,  we  found  that  38%  of 
the  execution  time  was  spent  doing  Fast  Fourier  Transforms,  implemented  in  FOR¬ 
TRAN.  This  indicates  that  the  code  can  be  run  faster  using  array  hardware.  In  addi¬ 
tion,  this  also  indicates  that  the  technique  may  be  faster  than  transformation  to  the 
k-f  domain  for  excising  the  surface  wave  and  inverse  transformation  to  the  x-t 
domain. 

Since  Saatpilar  and  Canitez  (1988)  have  shows  that  phase  match  filtering  to 
remove  ground  roll  significantly  improves  stacked  sections,  our  focus  was  on  making 
the  technique  easier  to  use  and  understand.  Our  implementation  of  adaptive  phase 
matched  filter  technique  is  very  robust  since  only  simple  control  parameters  are 
required  and  since  it  is  well  behaved,  even  for  data  with  no  noticeable  dispersion.  With 
future  24-bit  field  data  acquisition,  the  subtractive  nature  of  the  coherent  noise  removal 
may  significantly  improve  marginally  collected  data,  without  distorting  the  underlying 
exploration  signal. 
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Table  1.  Matched  filter  control  parameters 

PASS  /o  UQ  fx  Ux 

_ (Hz)  (m/sec)  (Hz)  (m/sec) 

1  6  914  20  305 

2  8  1219  40  610 

3  8  610  40  305 


Table  2.  Matched  filter  control  parameters 

PASS  /0  UQ  f !  Ux 

_ (Hz)  (m/sec)  (Hz)  (m/sec) 

1  2  366  30  183 

2  5  914  30  305 

3  5  366  20  61 

4  10  2743  40  914 


Table  3.  Matched  filter  control  parameters 

PASS  f0  U0  fx  Ux 

_ (Hz)  (tn/sec)  (Hz)  (m/sec) 

1  7  320  40  290 

2  7  1829  40  1219 
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Fig.  1.  Input  data  sets  from  the  Oklahoma  test  site,  a)  dynamite  source,  group  array 
DynGrp;  b)  dynamite  source  1  m  group  array  DynNoGrp;  c)  vibrator  source,  group 
array  VibGrp;  d)  vibrator  source,  1  m  group  array  VibNoGrp.  All  traces  are  AGC’d 
with  a  500  ms  window. 
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Fig.  2.  Data  from  Oklahoma  test  site  high  pass  filtered  with  a  comer  at  32  Hz.  a) 
DynGrp,  b)  DynNoGrp,  c)  VibGrp,  and  d)  VibNoGrp  data  sets.  All  traces  are 
AGC’d  with  a  500  ms  window. 
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Fig.  3.  Phase  matched  filter  processing  results  for  data  set  DynGrp.  a)  Original  data 
set,  b)  output  of  pass  1,  c)  output  of  pass  2,  d)  output  of  pass  3.  All  traces  are  AGC’d 
with  a  500  ms  window. 


Fig.  4.  Phase  matched  filter  processing  results  for  data  set  VibNoGrp.  a)  Original 
data  set,  b)  output  of  pass  1,  c)  output  of  pass  2,  d)  output  of  pass  3.  All  traces  are 
AGC’d  with  a  500  ms  window. 
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DynGrp,  b)  DynNoGrp,  c)  VibGrp,  and  d)  VibNoGrp  data  sets.  All  traces  are 
AGC’d  with  a  500  ms  window. 
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Fig.  6.  Results  of  passing  a  high  pass  filter  with  a  comer  at  32  Hz  to  the  phase 
matched  filter  output  of  Figure  5.  a)  DynGrp,  b)  DynNoGrp,  c)  VibGrp,  and  d)  Vib- 
NoGrp  data  sets.  All  traces  are  AGC’d  with  a  500  ms  window. 
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Fig.  7.  a)  Initial  west  Texas  record  section,  b)  result  of  a  4  passes  of  phase  match 
filters.  All  traces  are  AGC’d  with  a  500  ms  window. 
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Fig.  8.  a)  Initial  permafrost  record  section,  b)  result  of  2  passes  of  phase  match  filters. 
All  traces  are  AGC’d  with  a  500  ms  window. 
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ABSTRACT 


Source  parameters  of  59  large  intra-continental  earthquakes  as  given  in  recent  litera¬ 
ture  are  used  to  derive  seismic  scaling  relations  using  constrained  least-squares  for  the 
purpose  of  comparison  with  previously  proposed  scaling  relations  for  eastern  North  Amer¬ 
ica.  The  world-wide  data  of  intra-continental  earthquakes  are  consistent  with  a  power  3.0 
spectral  scaling  (constant  stress  drop)  as  previously  suggested  by  Somerville  et  al.  (1987) 
from  a  much  smaller  data  set  On  the  other  hand,  the  data  do  not  preclude  a  power  3.3 
spectral  scaling.  Consequently,  this  interpretation  permits  stress  drop  to  increase  with 
seismic  moment,  but  with  a  smaller  slope  and  at  a  lower  level  than  suggested  by  the 
revised  mid-plate  scaling  relation  of  Nuttli  et  al.  (1989).  Both  inteipretations  of  the  data 
predict  very  similar  comer  periods.  Using  the  proposed  scaling  relations,  magnitudes 
determined  from  random  process  theory  and  finite  source  modeling  compare  reasonably 
well  with  observations  in  eastern  North  America.  In  addition,  the  estimated  source  param¬ 
eters  of  selected  large  historical  earthquakes  are  only  slightly  different  from  those  previ¬ 
ously  suggested. 


INTRODUCTION 

Eastern  North  America  has  experienced  significant  historical  earthquakes  such  as  that 
of  1663  in  the  St.  Lawrence  Valley,  that  of  1755  near  Cape  Ann,  Mass.,  those  of  1811- 
1812  near  New  Madrid,  Mo.,  and  that  of  1886  near  Charleston,  S.C.  (Coffman  and  von 
Hake.  1973).  No  major  damaging  earthquake,  however,  has  occurred  there  since  the  instal¬ 
lation  of  regional  networks.  Large  intra-continental  earthquakes,  which  are  several  hun¬ 
dred  kilometers  distant  from  any  plate  boundary  (Scholz  et  al. ,  1986),  occur  infrequently, 
but  have  a  high  potential  for  destruction.  In  order  to  estimate  ground  motions  due  to  such 
large  events,  seismic  scaling  relations  must  be  employed.  These  relations  enable  an  esti¬ 
mate  of  magnitude  or  seismic  moment  from  observed  fault  dimensions  (geologic  observa¬ 
tions  and  microseismicity  studies,  e.g.,  Wyss,  1979)  and  a  determination  of  peak  ground 
motion  parameters  (e.g.,  Boore,  1983;  Atkinson,  1984;  Joyner,  1984;  Boore  and  Atkinson, 
1987;  EPRI,  1988;  Herrmann  and  Jost,  1988;  Heaton  and  Hartzell,  1989).  The  basic 
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assumption  of  seismic  scaling  laws  is  that  the  physics  of  material  failure  is  independent  of 
the  size  of  an  event,  i.e„  lacking  a  scale  length  (Mandelbrot,  1983).  This  seems  to  be 
valid,  however,  only  for  certain  subsets  of  earthquakes,  i.e.,  shallow  plate  margin,  mid¬ 
plate,  or  deep  earthquakes. 

For  mid-plate  earthquakes,  several  seismic  scaling  relations  have  been  proposed  that 
were  based  on  data  from  eastern  North  America  (Nuttli,  1983;  Hasegawa,  1983;  Nuttli 
et  al. ,  1989).  However,  it  is  currently  a  matter  of  debate  whether  or  not  scaling  for  mid¬ 
plate  earthquakes  differ  from  that  for  plate  margin  earthquakes  (Haar  et  al. ,  1986;  Somer¬ 
ville  et  al. ,  1987;  Nuttli  et  al. ,  1989).  Due  to  the  lack  of  large  instrumentally  recorded 
events,  a  decision  on  this  issue  for  North  America  is  hard  to  reach,  thus  making  it  difficult 
to  evaluate  seismic  hazard  (Boore  and  Atkinson,  1987;  Somerville  et  al. ,  1987;  Cop¬ 
persmith  and  Youngs,  1989). 

In  general,  four  essential  differences  are  noted  between  large  earthquakes  in  eastern 
North  America  and  those  of  western  North  America  (e.g.,  San  Francisco,  1906).  First,  the 
felt  areas  in  eastern  North  America  are  fundamentally  larger  than  those  in  the  west.  Hein¬ 
rich  (1941)  pointed  out  that  the  Charleston,  Missouri  earthquake  of  1895  caused  only 
moderate  damage  in  the  epicentral  region,  but  was  felt  as  far  away  as  the  Atlantic  coast, 
the  Gulf  coast,  New  Mexico,  and  Canada.  Gutenberg  and  Richter  (1949)  attributed  the 
wider  area  of  perceptibility  for  given  magnitude  to  greater  focal  depth.  This  conclusion 
was  contradicted  by  recent  data  that  show  the  maximum  focal  depth  for  earthquakes  in  the 
central  United  States  is  about  25  km  (Nuttli  and  Herrmann,  1987).  Nuttli  (1973b)  pro¬ 
posed  another  explanation:  he  inferred  that  the  coefficient  of  anelastic  attenuation  for  1  -Hz 
Lg- waves  is  about  0.0006  km~x  for  eastern  North  America,  as  compared  to  0.005  km~x 
for  coastal  California.  Further  studies  (e.g.,  Johns  et  al. ,  1977;  Bollinger,  1979;  Singh  and 
Herrmann,  1983)  on  the  coefficient  of  anelastic  attenuation,  or  apparent  Q,  in  the  central 
and  eastern  United  States  confirmed  Nuttli ’s  results  which  completely  explain  the 
differences  in  areas  of  perceptibility.  By  using  this  attenuation  relation,  Nuttli  (1973b) 
was  able  to  present  a  reliable  method  of  Lg-wave  magnitude  determination  (Baker,  1967). 

Second,  none  of  the  historical  earthquakes  in  eastern  North  America  has  produced 
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observable  surface  rupture  (Nuttli  and  Herrmann,  1987).  Large  events  like  those  near 
Charleston  in  1886  or  near  New  Madrid  in  1811  -  1812  produced  numerous  sand-blows, 
and  the  latter  even  rapids  in  the  Mississippi  river  due  to  subsidence  (Street  and  Nuttli, 
1984);  however,  surface  traces  of  faulting  comparable  to  those  observed  in  the  western 
United  States  have  not  been  found. 

Third,  assuming  that  magnitude  estimates  of  historical  events  are  not  grossly  in 
error,  earthquakes  in  eastern  North  America  seem  to  have  significantly  smaller  fault 
dimensions  than  in  western  North  America  for  a  given  magnitude.  For  example,  standard 
scaling  relations  require  fault  lengths  between  300  to  1000  km  for  earthquakes  with  Ms  = 
8.0  to  8.8  (e.g.,  Kanamori  and  Anderson,  1975;  Nuttli,  1985). 

Finally,  recurrence  rates  are  low  in  eastern  North  America  compared  to  the  west 
(Bakun  et  al.,  1986;  Johnston  and  Nava,  1984,  1985).  In  general,  repeat  times  are  related 
to  tectonic  loading  and  the  strength  of  preexisting  zones  of  weakness  (Talwani,  1989). 
Since  the  historical  record  is  significantly  smaller  than  one  complete  recurrence  cycle  of 
large  earthquakes,  future  large  events  may  occur  at  somewhat  unexpected  locations,  e.g., 
Meers  fault,  Oklahoma  (Hinze,  1988;  Coppersmith  and  Youngs,  1989;  Mitchell  et  al. , 
1989). 

The  last  two  features  of  large  earthquakes  in  eastern  North  America  imply  that  long 
repeat  times  correspond  to  short  fault  rupture  lengths.  Kanamori  and  Allen  (1986),  extra¬ 
polating  laboratory  experiments,  proposed  that  frictional  strength  on  the  fault  plane 
increases  as  the  contact-time  increases,  which  explains  the  observed  higher  stress  drops  for 
intraplate  events  (Kanamori  and  Anderson,  1975).  Put  differently,  events  with  long  repeat 
times  seem  to  have  a  larger  number  of  asperities  than  events  with  short  repeat  times 
(Kanamori  and  Allen,  1986).  The  same  conclusion  was  reached  by  Scholz  et  al.  (1986) 
who  inferred  that  large  intraplate  earthquakes  have  about  6  times  larger  slip  than  interplate 
events  of  the  same  size,  indicating  a  6  times  larger  stress  drop.  Previously,  Madariaga 
(1979)  indicated  that  larger  stress  drops  for  a  given  seismic  moment  indicate  a  rougher 
fault  in  agreement  with  the  above  arguments. 

Nuttli  (1983)  concluded  that  mid-plate  earthquakes  in  eastern  North  America  arc 


characterized  by  source  spectra  for  which  the  comer  period,  T 02,  is  proportional  to  the 
one-fourth  power  of  the  seismic  moment,  M0,  for  earthquakes  with  seismic  moments 
between  1020  and  1028  dyne-cm  (see  also  Street  and  Turcotte,  1977;  Iio,  1986).  This  scal¬ 
ing  requires  an  increasing  stress  drop  for  increasing  seismic  moment.  Haar  et  al.  (1986) 
questioned  Nuttli’s  (1983)  conclusion.  They  inferred  that  the  comer  period  is  proportional 
to  the  one-third  power  of  the  seismic  moment  for  both  western  and  eastern  North  Ameri¬ 
can  earthquakes  of  seismic  moments  between  1017  and  1026  dyne-cm.  This  scaling 
renders  the  stress  drop  constant,  independent  of  seismic  moment  (Kanamori  and  Anderson, 
1975).  Constant  stress  drop  scaling  is  a  result  of  assuming  a  similarity  condition  between 
large  and  small  events  (Aki,  1967),  i.e.,  fault  area  proportional  to  (fault  length)2,  u  pro¬ 
portional  to  fault  length  (constant  strain  drop;  Kanamori  and  Anderson,  1975;  Scholz, 
1982),  and  rupture  velocity  independent  of  seismic  moment. 

Somerville  et  al.  (1987)  related  source  duration  to  seismic  moment  for  earthquakes 
in  eastern  and  western  North  America  and  for  earthquakes  in  continental  interiors.  They 
inferred  a  logarithmic  slope  of  three  for  all  three  classes  of  earthquakes.  In  addition,  they 
found  that  the  spectral  scaling  relations  of  Nuttli  (1983)  predicted  source  durations  sub¬ 
stantially  larger  than  they  obtained  from  waveforms  of  short-period  P  waves.  Part  of  this 
incompatibility  may  be  related  to  differences  of  duration  of  P-waves  versus  S-waves  or 
higher  mode  surface  waves  (Hanks  and  Wyss,  1972;  Savage,  1972;  Cohn  et  al. ,  1982; 
Hasegawa,  1983;  Nuttli,  1983;  Silver,  1983).  Scholz  et  al.  (1986)  compared  large  inter- 
plate  and  intraplate  earthquakes  and  concluded  that  though  stress  drops  are  different,  both 
populations  follow  the  same  scaling  relation.  In  addition,  recent  data  suggest  that  relating 
a  comer  frequency  of  1.0  Hz  to  a  seismic  moment  of  1 022  dyne-cm  is  inappropriate  (e.g.. 
Shin  and  Herrmann,  1987;  Chun  et  al. ,  1989). 

The  primary  purpose  of  this  study  is  to  develop  scaling  relations  for  medium  to  large 
intra-continental  earthquakes  based  on  a  world-wide  data  set  and  to  compare  these  rela¬ 
tions  to  the  proposed  scaling  laws  for  eastern  North  America  (c.f.,  forthcoming  EPRI 
study  on  maximum  earthquakes.  Coppersmith  and  Youngs,  1989).  In  the  following,  we 
first  analyze  the  revised  mid-plate  scaling  law  (Nuttli  et  al. ,  1989).  Second,  recent  litera¬ 
ture  is  searched  for  source  parameters  of  large  intra-continental  earthquakes  world-wide. 
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These  source  parameters  are  used  to  derive  seismic  scaling  relations  which  are  constrained 
to  be  consistent  with  fundamental  seismic  laws.  Finally,  these  relations  are  compared  to 
the  revised  mid-plate  scaling  relation  (Nuttli  et  al. ,  1989),  the  scaling  relation  for  eastern 
Canada  (Hasegawa,  1983),  and  the  scaling  law  of  Somerville  et  al.  (1987).  Applying  the 
intra-continental  scaling  relations  for  eastern  North  America,  estimates  of  source  parame¬ 
ters  of  selected  large  historical  events  are  discussed. 

REVISED  MID-PLATE  SCALING  RELATION 

In  the  following,  we  give  an  analysis  of  the  recently  revised  mid-plate  scaling  na¬ 
tion  (Nuttli  et  al. ,  1989).  The  reduced  far-field  displacement  spectra  were  constructed  by 
trial  and  error  to  fit  the  relation  between  comer  frequency  and  seismic  moment,  observa¬ 
tional  data  of  body-wave  magnitude  versus  surface-wave  magnitude,  and  of  body-wave 
magnitude  and  surface-wave  magnitude  versus  seismic  moment  In  addition,  the  scaling 
relations  were  constrained  to  give  reasonable  values  of  the  fault  length  and  width,  of  aver¬ 
age  fault  displacement,  and  of  static  stress  drop  versus  seismic  moment.  In  general,  the 
determination  of  source  spectra  trades  off  with  frequency  dependent  anelastic  attenuation, 
site  amplifications,  instrument  influence  (band-limitation),  and  radiation  pattern  effects 
(Chun  et  al. ,  1989). 

Figure  1  displays  the  revised,  reduced  far-field  displacement  spectra  by  Nuttli  et  al. 
(1989)  which  show  the  following  characteristics:  At  long  periods,  each  spectmm  parallels 
the  period  axis  (slope  =  0);  the  level  is  proportional  to  the  seismic  moment  (Keilis-Borok, 
1959).  At  intermediate  periods,  spectra  of  small  earthquakes  have  one  comer  period  T02, 
where  the  slope  of  the  spectrum  changes  from  0  to  2  (assuming  the  co2  model,  Aki,  1967). 
This  comer  period  is  proportional  to  the  linear  dimension  of  the  fault  (Savage,  1972). 
Larger  earthquakes  can  have  two  comer  periods  Tu  and  T0 1(  the  shorter  proportional  to 
the  width,  the  longer  to  the  length  of  the  fault  plane,  respectively.  The  slope  of  the  spec¬ 
trum  between  the  two  comer  periods  is  unity,  indicating  a  non-uniform  stress  release  on 
the  fault  plane  or  a  partial  stress  drop  (Brune  1970,  1971;  Savage,  1972;  Brune  et  al. , 
1986).  The  short-period  trend  of  the  spectra  follows  an  (o-power  law  (Aki,  1967;  Brune, 
1970,  1971).  Data  from  eastern  North  America  seem  to  follow  the  ©-square  model  (Aki, 
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1967;  Chun  et  al. ,  1989),  and  not  the  o>-cube  model  (Aki,  1967;  Savage,  1972). 

The  core  of  a  scaling  relation  is  a  function  describing  the  increase  of  comer  period 
(for  larger  earthquakes  a  comer  period  7'02  can  be  simply  constructed)  with  seismic 
moment  (similarity  condition).  By  accepting  this  relation,  which  is  determined  from  small 
to  medium  size  events  in  eastern  North  America,  for  all  earthquakes,  it  is  possible  to 
predict  reduced  far-field  displacement  spectra  for  earthquakes  of  any  given  size.  The 
revised  mid-plate  scaling  law  is  described  by  the  following  relations  between  comer 
period  and  seismic  moment  (Figure  2).  For  fog  M0  <  23  : 


log  T Q2  =  -5.872  +  0.251  log  M0  ,  (la) 

or  for  the  comer  frequency  fc  =  I/T02 

fc  =  7.45  x  105  Mq 1/4  ,  (lb) 

where  comer  periods  are  in  sec,  comer  frequencies  in  Hz,  and  seismic  moments  in  dyne- 

cm.  Note  that  T01  =  T02  =  Tn  for  fog  M0  <  23  rendering  the  aspect  ratio  (fault 

width/fault  length)  constant  for  small  events.  For  larger  seismic  moments,  the  revised 

mid-plate  scaling  relation  causes  the  aspect  ratio  to  differ  from  unity  (Kanamori  and 


Anderson,  1975;  Heaton  and  Hartzell,  1988). 

log  T01  =  -7.214  +  0.309  log  M0 
log  T 02  =  -6.373  +  0.273  log  M0 
or  for  the  comer  frequency  fc 

fc  =  2.36  x  106  Afo 1/3  67 
log  Tn  =  -5.455  +  0.233  log  M0 
For  log  M0  >  24  : 

log  Tq,  =  -8.033  +  0.343  log  M0 
log  Tm  =  -6.589  +  0.282  log  M0 
or  for  the  comer  frequency  fc 

fc  =  3.88  x  106  A/o1/3'55 
log  T, 2  =  -5.155  +  0.220  log  M0 
The  relation  between  comer  frequency  and 

(1983)  formula  (20  <  log  M0  S  28): 

fc  =  3.55  x  105  MqVA  , 


For  23  £  fog  M  0  <  24  : 

(2) 

(3a) 


(3b) 

(4) 


(5) 

(6a) 


(6b) 
(7) 

seismic  moment  can  be  compared  to  Nuttli’s 


(8) 
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and  that  for  plate  margin  events  assuming  a  constant  stress  drop  of  100  bar  (Boore  and 
Atkinson,  1987) 

fc  =  7.96  x  107  M  o1/3  .  (9) 

Furthermore,  it  can  be  compared  to  the  relation  by  Atkinson  (1984,  equation  (11))  which 

is  based  on  large  eastern  Canadian  events  (23.9  <  log  Mq  <  25  8): 

fc  =  2.00  x  106  Mq1/3-51  .  (10) 

Note  the  close  agreement  with  the  revised  mid-plate  scaling  law  (6b). 

The  revised  mid-plate  scaling  relation  can  also  be  described  by  (Figure  1): 


log  M0  =  1.0  mb  +  17.5  for  mb  <  4.2  (1 1) 

log  Mq  =  2.0  mb  +  13.1  for  mb  >  4.2  (12) 

log  M0  =  1.0  Ms  +  18.85  for  Ms  <  8.15  (13) 

log  Mq  =  4.0  log  +  23.4  for  <  Is  (14) 

log  Mq  =  3.55  log  T 02  +  23.4  for  >  1  s.  (15) 


Table  1  gives  values  of  body-wave  magnitude,  surface-wave  magnitude,  moment 
magnitude,  and  the  comer  periods  Tox,  T 12  and  T0 2  for  seismic  moments  between  1020 
and  1028  dyne-cm.  It  also  contains  estimated  values  of  the  fault  rupture  length  (measured 
along  strike)  and  width  (measured  along  dip),  the  static  stress  drop,  and  the  average  fault 
displacement  that  are  calculated  from  equations  given  in  Brune  (1970,  1971),  Savage 
(1972),  and  Geller  (1976): 


L  =  3.57  p  T0ll 2k 

(16) 

W  =  2.30  v*  r12/7t 

(17) 

Ao  =  7  Mq/  16  (LW/k)3'2 

(18) 

Mq  =  p  L  W  u 

(19) 

In  the  following,  fault  length  and  width,  L  and  W,  are  in  cm;  the  comer  periods  in 
seconds;  the  shear-wave  velocity,  (5,  and  the  rupture  velocity,  vR ,  in  cm/s ;  the  stress  drop 
Ao  in  bar,  the  average  displacement,  u,  in  cm;  and  the  shear  modulus,  p,  in  dyne/cm2. 
The  rupture  velocity  is  assumed  to  be  constant  and  equal  to  0.9  x  p.  The  value  of  p  is 
assumed  to  be  3.3  x  1011  dynes/cm2.  Bonilla  et  al.  (1984)  estimated  that  p  can  vary 
between  1.7  x  10n  dyne/cm2  and  3.4  x  1011  dynelcm2. 

Following  Kanamori  (1978)  and  Hanks  and  Kanamob  (1979),  the  moment  magni- 


tude  Mw  is  given  by 

Mw  =  j  log  M0  -  10.73  ,  (20) 

which  was  introduced  to  extrapolate  the  surface-wave  magnitude  scale  beyond  its  satura¬ 
tion  point  at  about  Ms  =  8.0  (e.g.,  Figure  7  in  Geller,  1976;  Nuttli,  1985).  On  the  other 
hand,  Singh  and  Havskov  (1980)  pointed  out  that  the  equivalence  between  values  of  Mw 
and  Ms  for  Ms  <  8.0  is  valid  only  for  interplate  events.  For  intraplate  earthquakes,  Singh 
and  Havskov  (1980)  proposed  a  correction  term  of  0.27  to  the  moment  magnitude  as 
given  in  (20)  rendering 

Mw  =  —  log  Mq  -  10.46  .  (21) 

In  Table  1,  the  definition  (21)  is  used.  However,  the  correspondence  between  Ms  and  Mw 
is  poor  for  the  revised  mid-plate  scaling  relation  (Table  1).  This  fit  is  not  improved  if  (20) 
instead  of  (21)  is  used. 

Table  2  shows  similarity  conditions  (Kanamori  and  Anderson,  1975)  and  rise  times 
(Geller,  1976)  for  the  revised  mid-plate  scaling  relation.  For  small  events,  the  aspect  ratio 
is  one  and  decreases  for  larger  events  as  the  fault  length  increases  with  respect  to  fault 
width.  The  strain  drop  («/ L)  increases  with  increasing  moment  indicating  that  large  mid¬ 
plate  events  release  more  strain  than  small  events.  Finally,  the  dynamic  similarity 
decreases  slightly,  implying  somewhat  different  effective  stress  for  small  and  large  events. 

Nuttli  et  al.  (1989)  pointed  out  that  the  revised  scaling  relation  (Table  1)  gives  a 
slightly  better  fit  to  observational  data  of  mb  versus  Ms,  mb  versus  M0,  and  Ms  versus 
Mq  than  the  earlier  spectral  scaling  model  of  Nuttli  (1983).  The  improvement  is  significant 
for  earthquakes  of  moment  less  than  1 022  dyne-cm.  For  larger  earthquakes,  the  predicted 
curves  given  by  Figure  1  and  Table  1  are  almost  identical  to  the  original  (Nuttli,  1983) 
curves.  The  fit  of  the  comer  period  (T02)  data  (e.g.,  Fletcher,  1982;  Shin  and  Herrmann, 
1987;  Somerville  et  al. ,  1987)  to  those  predicted  by  the  revised  spectral  scaling  relation 
(Figure  3)  seems  to  be  superior  to  the  fit  of  the  relations  by  Nuttli  (1983)  and  Somerville 
et  al.  (1987).  Hasegawa’s  (1983)  relation  is  essentially  identical  to  that  of  Nuttli  et  al. 
(1989).  Furthermore,  stress  drops  can  be  compared  to  source  duration  (Figure  4).  Again, 
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the  revised  mid-plate  scaling  relation  seems  to  produce  a  superior  fit  to  the  data  of  Somer¬ 
ville  et  al.  (1987)  and  of  Shin  and  Herrmann  (1987)  after  neglecting  stress  drops  larger 
than  200  bar.  However,  the  predicted  stress  drop  values  are  somewhat  larger  than  the 
values  given  in  Street  and  Turcotte  (1977,  Figure  4)  which  may  be  due  to  differences  in 
comer  period  estimates  (e.g.,  Haar  et  al. ,  1986;  Anderson,  1986;  Andrews,  1986).  From 
(18),  a  small  error  in  measuring  comer  periods  from  data  will  significantly  alter  the  stress 
drop  value  (Kanamori  and  Anderson,  1975;  Brune  et  al.,  1979;  Atkinson,  1984).  Note, 
that  the  presence  of  asperities  on  the  fault  plane  can  also  alter  the  comer  period  (e.g..  Fig¬ 
ure  6.16  in  Kostrov  and  Das,  1988). 

Recently,  Chun  et  al.  (1989)  investigated  source  spectra  of  Miramichi  earthquakes 
with  mig  =  2.6  -  5.4.  They  found  strong  support  for  an  increase  of  stress  drop  with 
seismic  moment  (Nuttli,  1983)  after  correlating  a  seismic  moment  of  5.6  x  1023  dyne-cm 
to  a  comer  frequency  of  1  Hz  following  Shin  and  Herrmann  (1987). 

On  the  other  hand.  Hanks  (1982)  and  Boore  (1986b)  pointed  out  that  the  observed 
increase  of  stress  drop  with  seismic  moment  (breakdown  in  similarity)  for  earthquakes 
with  magnitudes  less  than  about  5  could  be  due  to  a  process  that  limits  high  frequencies 
(f  mu)>  where  it  is  irrelevant  whether  the  process  is  due  to  instrument,  source,  or  site. 
This  band-limiting  effect  can  lead  to  significant  differences  in  scaling  for  large  and  small 
earthquakes  (Boore,  1986b).  If  such  an  effect  is  present  in  eastern  North  America  where 
there  are  no  data  for  large  events,  the  extrapolation  of  properties  of  small  to  large  events 
can  be  in  error.  It  is  interesting  to  note  that  the  revised  mid-plate  scaling  relation  (Nuttli 
et  al. ,  1989)  indeed  shows  different  scaling  for  small  and  large  earthquakes  (Figure  2). 

SCALING  RELATIONS  FOR  INTRA- CONTINENTAL  EARTHQUAKES 

To  test  the  revised  mid-plate  scaling  relation  derived  from  earthquakes  in  eastern 
North  America,  a  search  of  the  recent  literature  for  observed  source  properties  of  world¬ 
wide  intra-continental  earthquakes  was  performed.  In  this  study,  we  restrict  intraplate 
earthquakes  as  recently  reviewed  by  Talwani  (1989)  to  continental  areas.  However,  intra¬ 
continental  events  are  understood  to  be  less  restrictive  than  events  in  stable  continental 
interiors  as  defined  by  Johnston  (1989).  In  this  way,  we  could  include  high  quality  data 
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for  events  in  central  Asia  (Figure  5)  which  is  considered  an  active  plate  interior  by  Johns¬ 
ton  (1989,  his  Figure  1).  On  the  other  hand,  none  of  the  events  chosen  were  located  on 
the  Himalayan  Frontal  Thrust,  the  plate  boundary  between  India  and  Asia.  For  example, 
we  omitted  the  events  Bihar-Nepal  in  1934  and  Assam  in  1950  (Chen  and  Molnar,  1977). 
We  also  omitted  events  in  the  Garm  region  (e.g.,  Markansu,  1974;  Jackson  et  al. ,  1979; 
Hamburger  et  al. ,  1988).  Table  3  gives  the  hypocenters  and  fault  plane  solutions  of  the 
events  used  and  Table  4  the  source  parameters.  In  both  tables,  the  source  of  the  data  is 
indicated.  All  events  chosen  were  located  in  the  crust. 

In  Table  4,  fault  lengths  determined  from  observed  fault  scarps  or  fault  lengths 
determined  by  waveform  matching  (e.g.,  Cipar,  1979)  were  preferred  over  values  deter¬ 
mined  from  aftershock  studies.  Darragh  and  Bolt  (1987)  pointed  out  that  there  can  be  a 
discrepancy  between  the  extent  of  the  aftershock  zone  and  the  surface  scarp  of  faulting. 
Zhou  et  al.  (1983b)  concluded  that  there  are  earthquakes  (e.g.,  1967  Ganzi,  1975 
Haicheng)  where  the  aftershock  zone  significantly  overestimates  the  fault  length.  Cipar 
(1979)  and  Shedlock  et  al.  (1985)  found  the  aftershock  zone  for  the  Haicheng  earthquake 
significantly  larger  than  the  surface  trace  of  faulting.  Furthermore,  aftershock  zones  vary 
with  time  (Upadhyay  and  Duda,  1980)  and  their  location  might  not  be  well  constrained. 
On  the  other  hand,  surface  rupture  might  not  be  observable  (e.g.,  events  in  eastern  North 
America)  or  if  observable  (e.g.,  in  Australia)  might  not  extend  over  the  whole  length  that 
actually  ruptured  below  surface  even  for  shallow  events  (Wyss,  1979;  Bonilla  et  al. , 
1984).  Bonilla  et  al.  (1984)  argued  that  the  surface  rupture  length  can  be  taken  as  fault 
length  for  steeply  dipping  faults  and  aspect  ratios  smaller  than  0.5. 

Most  of  the  average  displacement  values  of  Table  4  are  from  field  observations. 
Average  slips  measured  at  the  earth’s  surface  show  a  large  degree  of  variability  especially 
for  the  very  largest  events  (e.g.,  Deng  et  al. ,  1986;  Zhang  et  al. ,  1987).  This  makes  the 
determination  of  an  average  value  of  displacement  along  the  surface  fault  trace  difficult.  In 
addition,  it  is  then  assumed  that  the  observed  displacements  at  the  earth’s  surface  are 
representative  for  the  whole  fault  plane. 

Rupture  widths  as  given  in  Table  4  are  determined  from  the  distribution  of  aft- 
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ershocks  or  microseismicity  studies  indicating  the  extent  of  the  seismogenic  zone. 
Fredrich  et  al.  (1988)  argued  for  taking  the  lower  bound  of  Australian  earthquakes  at  2 
times  the  centroid  depth. 

Most  duration  values  given  in  Table  4  were  determined  from  waveform  matching 
and  seem  to  follow  the  definition  by  Cohn  et  al.  (1982). 

Linear  Regression  Analysis 

Following  Marie  (1977),  a  linear  regression  of  source  parameter  Y  (e.g.,  logarithm  of 
fault  length)  on  log  M0  is  written  as  Y  =  A  +  B  log  M0.  This  equation  passes  through  the 
most  probable  value  of  Y  for  any  given  log  M0  and  this  equation  is  to  be  used  for 
estimating  Y  from  a'giyen  seismic  moment. 

To  obtain  the  most  probable  seismic  moment  from  a  given  source  parameter  Y,  we 
need  to  do  the  regression  log  M0  =  A'  +  B'  Y.  The  use  of  this  equation  for  estimating  Y 
from  a  given  seismic  moment  can  lead  to  noticeable  errors  (Mark,  1977;  Bonilla  et  al. , 
1984). 

In  this  study,  we  accept  the  seismic  moment  as  the  fundamental  parameter  of  earth¬ 
quake  size,  consequently  using  the  former  regression  analysis.  The  data  are  least  squares 
fit  to  a  straight  line  (Singh  et  al. ,  1980;  Bonilla  et  al. ,  1984)  using  a  routine  by  Press 
et  al.  (1986).  In  that  routine,  the  independent  variable  (log  M0)  is  assumed  error  free 
(Bolt,  1978)  and  a  standard  deviation  is  associated  with  Y  acting  as  a  weight.  We  feel  we 
are  not  able  to  assign  objective  weights  to  the  observations  in  Table  4  since  the  number  of 
independent  observations  of  source  parameters  for  each  earthquake  is  not  significantly 
larger  than  1  (Coppersmith  and  Youngs,  1989).  Therefore,  we  assume  a  constant  standard 
deviation  of  0.3  in  the  logarithmic  domain  for  all  points,  which  is  a  conservative  estimate 
for  most  source  parameters  (Wyss,  1979;  Bonilla  et  al. ,  1984). 

Magnitude  Relations 

Body-wave  magnitudes  from  Table  4  are  shown  in  Figure  6.  The  observed  scatter 
leads  to  question  the  usefulness  of  the  mb  magnitude  scale  for  the  classification  of  earth¬ 
quakes  according  to  size.  One  explanation  of  the  observed  scatter  in  Figure  6  might  be 


that  observed  mb  values  are  not  strongly  related  to  the  size  (i.e.,  seismic  moment)  of  an 
earthquake  (Dziewonski  and  Woodhouse,  1983).  Since  both  estimates  of  the  size  of  an 
earthquake  are  from  different  parts  of  the  radiated  spectrum,  these  differences  might  indi¬ 
cate  spectral  variations  of  source  radiation  as  suggested  by  Duda  and  Nortmann  (1983). 
On  the  other  hand,  Hanks  (1979)  showed  that  differences  in  radiated  spectral  amplitudes 
need  not  translate  into  differences  of  the  1  Hz  time  domain  amplitudes  used  for  determin¬ 
ing  mb  due  to  the  characteristics  of  energy  release  during  the  duration  of  faulting. 
Another  cause  of  the  scatter  seen  in  Figure  6  might  be  related  to  the  rather  inhomogene¬ 
ous  sources  of  mb  values  (differences  in  observatory  practice;  e.g.,  Willmore,  1979;  Aki 
and  Richards,  1980  (Appendix  2);  Houston  and  Kanamori,  1986;  Boore,  1986a;  McLaugh¬ 
lin  and  Jih,  1988).  In  addition,  mb  values  of  dip-slip  earthquakes  are  approximately  0.3 
magnitude  units  larger  than  for  strike-slip  earthquakes  (Houston  and  Kanamori,  1986).  In 
conclusion,  the  above  arguments  seem  to  suggest  that  the  mb  magnitude  scale  is  not  very 
useful  for  classifying  large  earthquakes. 

Accepting  seismic  moment  as  fundamental  parameter  of  earthquake  size,  we  per¬ 
formed  a  regression  of  observed  mb  magnitudes  on  log  M0,  assuming  a  standard  deviation 
of  0.3  for  each  magnitude  value. 

mb  =  0.23(±1.26)  +  0.22(±0.05)  log  M0  .  (22a) 

The  correlation  coefficient  is  0.63  indicating  a  poor  fit  of  the  data  to  a  straight  line.  This 

poor  correlation  leads  to  a  significantly  different  relation  if  we  take  mb  as  the  fundamental 

parameter  of  earthquake  size,  performing  a  regression  of  log  M0  on  mb  (written  in  reverse 

order) 

mb  =  -8.52  +  0.56  log  M 0  .  (22b) 

(22a)  is  significantly  different  from  the  relation  (12)  by  Nuttli  et  al.  (1989)  for  mb  >  4.2 
(short  dashes  in  Figure  6): 

mb  =  -6.55  +  0.50  log  M0  ,  (23) 

or  Hasegawa’s  (1983)  relation  for  4.2  <  M  <  6.6 

M  =  -7.30  +  0.54  log  Mq  ,  (24) 

where  M  =  mN  (i.e.,  mLg)  or  derived  from  ML.  Note  the  similarity  of  (22b),  (23),  and 
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(24). 


In  this  study,  we  estimate  the  most  likely  magnitude  for  a  given  seismic  moment, 
whereas  Nuttli  et  al.  (1989)  and  Hasegawa  (1983)  estimated  the  most  likely  seismic 
moment  from  a  given  magnitude.  From  Figure  6,  we  see  that  Nuttli  et  al.  (1989)  did  not 
consider  the  saturation  of  the  mb  scale  (e.g.,  Figure  1  in  Boore,  1986a).  On  the  other 
hand,  (22a)  seems  also  biased  due  to  including  the  saturation,  hence  predicting  too  large 
mb  values  at  smaller  seismic  moments.  A  change  of  slope  in  the  body-wave  magnitude 
relation  at  about  log  M0  =  24.5  seems  necessary  to  account  for  the  saturation  of  the  mb 
scale. 

The  surface-wave  magnitudes  of  Table  4  are  shown  in  Figure  7.  For  the  linear 
regression  of  Ms  on  log  M0,  magnitude  values  larger  than  8.0  and  moments  larger  than 
1028  dyne-cm  were  rejected,  and  a  standard  deviation  of  0.3  magnitude  unit  was  assigned 
to  each  magnitude. 

Ms  =  -13.87(±1.14)  +  0.79(±0.04)  log  M0  ,  (25) 

which  can  be  compared  to  (21)  (Singh  and  Havskov,  1980),  assuming  Mw  =  Ms  for  L  < 

100  km.  Comparing  magnitude  values  predicted  by  (25)  (solid  line  in  Figure  7)  and  (21) 
Gong  dashes  in  Figure  7)  in  the  moment  range  of  1025  and  1028  dyne-cm  shows  that  mag¬ 
nitude  differences  are  smaller  than  the  usual  standard  deviation  of  magnitude  observations 
of  0.3  magnitude  units.  On  the  other  hand,  (13)  of  the  revised  mid-plate  scaling  relation 
(Nuttli  et  al. ,  1989)  predicts  larger  Ms  values  for  the  same  moment  range  (short  dashes 
in  Figure  7).  If  we  reverse  the  order  of  the  linear  regression,  we  get  a  relation  very  simi¬ 
lar  to  (25),  predicting  magnitude  values  identical  to  those  from  (25)  in  the  moment  range 
studied.  Equally,  the  correlation  coefficient  is  0.96  indicating  a  good  fit  of  the  data  to  a 
straight  line. 

Fault  Areas  and  Average  Displacements 

Fault  areas,  A  =  L  W,  and  average  displacements  on  the  fault,  u,  taken  from  Table  4 
are  displayed  in  Figure  8  and  9,  respectively.  In  the  linear  regression  routine,  each  loga¬ 
rithmic  value  of  fault  area  or  average  displacement  was  assigned  a  standard  deviation  of 
0.3.  For  fault  areas,  a  factor  of  2  uncertainty  is  a  conservative  estimate  (Geller,  1976);  for 
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average  displacements,  it  may  actually  be  on  the  low  side.  To  be  consistent,  we  use  cm 
as  unit  of  length. 

log  A  =  -2.872(±1.544)  +  0.593(±0.058)  log  M0  (26) 

log  H  =  -6.730(±1.396)  +  0.337(±0.052)  log  M0  .  (27) 

If  we  neglect  the  standard  deviations,  these  separate  regressions  are  not  consistent  with 

(19).  Combined  and  constrained  regressions  can  be  used  to  force  agreement  with  (19). 
We  fit  the  data  of  fault  areas  and  average  displacements  simultaneously  to  two  different 
straight  lines,  imposing  the  following  constraints:  the  sum  of  the  slopes  had  to  be  unity, 

and  the  sum  of  the  intercepts  had  to  equal  -log  p  with  p  =  3.0  x  1011  dynes! cm1.  The 

corresponding  coefficients  had  somewhat  large  standard  deviations  indicating  that  the  data 
can  be  fit  by  a  multitude  of  relations.  Therefore,  we  decided  to  employ  further  constraints 
to  yield  two  sets  of  regression  relations.  For  the  first,  we  note  that  (26)  seems  compatible 
with  log  A  proportional  to  0.667  log  M0  of  Kanamori  and  Anderson  (1975).  Hence,  we 
also  constrain  the  slope  of  the  straight  line  for  the  area  relation  to  0.667.  This  first 
approach  leads  to  constant  stress  drop  scaling.  The  combined,  constrained  regression 
gives: 

log  A  =  -4.843  +  0.667  log  M0  (28) 

log  a  =  -6.631  +  0.333  log  M0  .  (29) 

For  the  second  approach,  we  note  that  fault  area  is  better  constrained  than  average 
fault  displacements,  e.g.,  the  correlation  coefficients  are  0.90  and  0.78  in  (26)  and  (27), 
respectively.  Hence  our  second  approach  assumes  a  slope  of  0.600  in  (26).  This  second 
approach  leads  to  variable  stress  drop  scaling.  The  combined,  constrained  regression 
gives 

log  A  =  -3.059  +  0.600  log  M0  (30) 

log  u  =  -8.415  +  0.400  log  M 0  .  (31) 

Note  that  the  results  of  the  combined,  constrained  regressions,  i.e.,  (28)  -  (31),  are 

modifications  of  (26)  and  (27)  that  are  within  the  range  of  their  standard  deviations. 

Stress  Drops 

From  (28),  an  equivalent  fault  radius  can  be  determined  assuming  Brune’s  model 
(1970,  1971).  Using  (18)  yields: 
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Ac  =  45  bar  (32) 

Hence,  the  data  in  Table  4  can  be  described  by  constant  stress  drop  scaling.  The  value  of 

the  stress  drop  is  somewhat  lower  than  suggested  by  Somerville  et  al.  (1987)  (Figure  10). 


Equivalently,  (30)  gives  with  (18) 

Ac  =  9.45  x  10"2  M001  ,  (33) 

indicating  increasing  stress  drop  (Figure  10)  as  proposed  by  Nuttli  et  al.  (1989).  How¬ 
ever,  (33)  gives  a  stress  drop  that  increases  slower  and  at  a  lower  level  than  suggested  by 
Nuttli  et  al.  (1989)  (Figure  10),  i.e.,  from  20  bar  at  1023  dyne-cm  to  60  bar  at  1028 
dyne-cm. 

The  observed  values  of  stress  drop  show  a  rather  large  scatter  of  an  approximately 
circular  shape  (the  correlation  coefficient  is  less  than  0.3)  and  are  somewhat  lower  than 
those  of  the  intra-plate  events  studied  by  Kanamori  and  Anderson  (1975).  However, 
Hasegawa  (1983)  noted  the  stress  drop  to  range  between  10  -  50  bar  for  4.2  <  M  <  6.6 
for  eastern  Canadian  events.  Both  relations,  (32)  and  (33),  seem  to  pass  closer  to  the 
center  of  that  circular  distribution  of  c  served  stress  drop  values  than  the  relations  by 
Somerville  et  al.  (1987)  or  Nuttli  et  al.  (1989). 

Durations 

Most  of  the  durations  given  in  Table  4  are  the  lengths  of  the  source  time  functions 
(triangular  or  trapezoidal  pulses)  as  defined  by  Cohn  et  al.  (1982).  We  excluded  very 
large  events  with  durations  longer  than  20  seconds  from  the  least  squares  fit  which  gives 
(Figure  11) 

log  x  =  -7.673(±1.547)  +  0.326(10.061)  log  M0  (34) 

with  a  correlation  coefficient  of  0.85.  From  Brune  (1970,  1971) 

log  M0  =  3  log  Tc  +  log  Ao  +  21.70  .  (35) 

For  constant  stress  drop  scaling,  we  substitute  (32)  and  (34)  into  (35)  assuming  x  =  Tc 

(Somerville  et  al. ,  1987).  This  consistency  check  gives  a  residual  of  0.33  in  the 
predicted  log  MQ  and  a  slight  difference  in  M0  dependence.  To  match  the  M0  depen¬ 
dence,  we  fix  the  slope  to  0.333  in  (34)  and  solve  in  the  regression  for  the  intercept, 
obtaining 
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log  T  =  -7.859  +  0.333  log  M0  ,  (36) 

which  gives  a  residual  of  -0.23,  that  can  be  completely  explained  by  a  shear  wave  velocity 


of  4.2  km/s  in  the  derivation  of  (35). 

Next,  we  investigate  increasing  stress  drop  scaling.  Substituting  (33)  and  (34)  into 
(35)  gives  a  residual  of  -2.34  and  a  slight  difference  in  M0  dependence.  To  match  the  M  0 
dependence  in  (35),  we  suggest 


log  t  =  -7.010  +  0.300  log  Mq  ,  (37) 

resulting  in  a  residual  of  -0.35. 

Iio  (1986)  and  Somerville  et  al.  (1987)  implied  that  duration  essentially  equals 
comer  period.  But  Herrmann  and  Goertz  (1981)  showed  that  a  symmetric  triangular  pulse 
of  duration  x  has  a  spectral  comer  period  of  n  x/2.  Assuming  that  all  durations  in  Table  4 
are  from  symmetric  triangular  pulses,  we  need  to  add  0.2  to  all  values  of  log  x  changing 
the  intercept  in  (37)  to  -6.814.  The  consistency  check  on  Brune’s  model  gives  then  a  resi¬ 
dual  of  +0.23,  which  could  be  completely  explained  by  a  shear  wave  velocity  of  3.0  km/s, 
which  is  rather  low. 


Another  way  of  testing  (36)  and  (37)  is  to  compare  it  to  predictions  of  the  duration 
as  given  by  Cohn  et  al.  (1982): 


R  ,  Rs in5  .  1 6R 
x  =  _  --  +  — - —  + 


0.8p  '  P  r  7icp  ’  (38) 

where  5  is  the  angle  between  the  normal  to  the  fault  plane  and  the  ray  path  to  the  station. 

The  first  two  terms  describe  the  rupture  propagation  time;  the  third  term  denotes  the  rise 
time  (Geller,  1976).  Assuming  p  =  3.5  km/s  and  <sin  5>  =  0.64  (Cohn  et  al. ,  1982),  and 
using  (28)  and  (30)  to  determine  the  equivalent  fault  radius  gives  relations  for  x  that  com¬ 
pletely  agree  with  (36)  and  (37),  respectively.  This  is  expected  since  Cohn  et  al.  (1982) 
also  followed  Brune  (1970,  1971). 


In  conclusion,  the  revised  duration  relation  (36)  suggests  a  power  3.0  spectral  scaling 
relation  as  previously  inferred  by  Somerville  et  al.  (1987)  from  a  much  smaller  o  a  set. 
On  the  other  hand,  the  data  seem  also  compatible  with  a  power  3.3  spectral  scaling  as 
indicated  by  (37).  Hasegawa  (1983)  inferred  for  large  eastern  Canadian  earthquakes 


log  Tc  =  -6.50  +  0.28  log  M  o  , 


(39) 
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which  is  essentially  identical  to  (10)  (Atkinson,  1984)  and  (6a)  (Nuttli  et  al. ,  1989). 
Hence,  both  (36)  and  (37)  are  different  from  the  power  3.6  spectral  scaling  relation  pro¬ 
posed  for  eastern  North  America. 

Fault  Lengths  and  Rupture  Velocities 

Fault  lengths  from  Table  4  are  plotted  in  Figure  12.  Conservatively  assuming  a 
standard  deviation  in  the  logarithmic  domain  of  0.3  for  each  data  point,  we  get  from  the 
linear  regression  (Figure  12) 

log  L  =  —4.681(±1.080)  +  0.425(±0.041)  log  M0  .  (40) 

The  correlation  coefficient  is  0.93.  Note  that  (40)  supports  proportionality  of  fault  length 

and  average  displacement  if  (31)  for  variable  stress  drop  scaling  is  used. 

For  estimating  rupture  velocities,  we  make  the  assumption  that  the  duration  can  be 
interpreted  as  x  =  R/vr  or  x  =  2 R/vr  for  bilateral  and  unilateral  rupture  propagation, 
respectively  (Sato  and  Hirasawa,  1973;  Boatwright,  1980;  Iio,  1986).  Using  (28)  and 
(36),  we  get  vR  =  1.5  km/s  and  3.1  km/s  for  bilateral  and  unilateral  rupture  propagation, 
respectively  (assuming  constant  stress  drop  scaling).  Using  (30)  and  (37),  we  get  1.7  km/s 
and  3.4  km/s  for  bilateral  and  unilateral  rupture  propagation,  respectively.  Note  that  those 
values  are  consistent  with  observed  rupture  velocities  that  vary  between  1.0  km/s  -  4.8 
km/s  (Purcaru  and  Berckhemcr,  1982). 

Using  a  rupture  length  as  determined  from  slightly  modifying  (40)  (L  a  M o'4) 
would  lead  to  vR  a  Mq1  .  For  a  bilateral  rupture  propagation,  vR  a  L025  which  would  be 
similar  to  Iio’s  (1986)  suggestion  of  vR  a  L020  -  L0,38  (see  also  Purcaru  and  Berckhemer, 
1982). 

Spectral  scaling 

The  scaling  relations  derived  above  are  used  in  an  effort  to  construct  a  spectral  scal¬ 
ing  model.  We  follow  the  approach  of  Nuttli  (1983)  and  Nuttli  et  al.  (1989)  and  take  vR 
=  0.9  p.  Assuming  the  to-squared  model,  it  is  sufficient  to  specify  the  comer  periods  to 
completely  describe  the  spectral  scaling. 
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For  constant  stress  drop  scaling,  we  follow  Savage  (1972),  i.e.,  (16),  and  using  (40) 

log  70!=  -9.980  +  0.425  log  M0  .  (41) 

From  (17),  (28),  and  (40),  we  find 

log  r12= -5.525  +  0.242  log  A/0  .  (42) 

Furthermore,  we  take  (36)  for  7’02.  A  quick  test  shows  that  the  relation  for  T 12  seems 

incompatible  with  (36)  and  (41).  We  choose  to  adapt  T 12  and  take  T 12  =  0.62  s  for  log 

M0  =  23.  Using  the  geometry  of  the  reduced  far-field  displacement  spectra,  we  get  T J2  = 

19.3  s  for  log  Mq  =  28.0.  Then 

log  r12=  -7.058  +  0.298  log  M0  .  (43) 

Figure  13a  shows  that  the  spectral  scaling  relations  (solid  lines)  predict  longer  comer 

periods  than  those  from  Nuttli  et  al.  (1989,  dashed  lines).  At  smaller  moments,  however, 

both  spectral  scaling  relations  predict  similar  comer  periods.  Included  in  Figure  13a  are 

comer  periods  estimated  from  one  station  for  events  #  40  and  #  41  (Upadhyay  and  Duda, 

1980).  The  scaling  relations  for  intra-continental  earthquakes  using  constant  stress  drop 

scaling  are  summarized  in  Table  5,  where  its  last  four  columns  will  be  discussed  below. 

Figure  13b  shows  a  comparison  of  scaling  relations  derived  for  the  case  of  constant 
and  variable  stress  drop  ( log  7"  12  =  -5.416  +  0.233  log  M0 ).  Note  that  both  relations 
predict  comer  periods  that  are  very  similar  (Table  6). 

Table  7a  shows  the  similarity  conditions  (Kanamori  and  Anderson,  1975)  and  rise 
times  (Geller,  1976)  for  constant  stress  drop  scaling.  For  small  events,  the  aspect  ratio  is 
one  and  decreases  for  larger  events  as  the  fault  length  increases  with  respect  to  fault 
width.  Values  of  the  aspect  ratio  were  derived  from  fault  length  and  fault  width  values  in 
Table  5.  Following  Savage  (1972),  note  that  T0 j  is  proportional  to  L,  and  Tn  proportional 
to  W.  Using  (41)  and  (43)  gives 

^  =  —  =  8.4  x  102  M^0127  ,  (44) 

L  T'oi 

which  gives  similar  values  for  the  aspect  ratio. 

The  strain  drop  decreases  with  increasing  moment  (Table  7a).  The  dynamic  similar¬ 
ity  decreases  slightly,  similar  to  the  relation  by  Nuttli  et  al.  (1989),  implying  somewhat 
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different  effective  stress  for  small  and  large  events.  Table  7b  shows  similarity  conditions 
for  the  case  of  variable  stress  drop.  Note  that  the  strain  drop  remains  constant. 


DISCUSSION  AND  APPLICATION  TO 
EASTERN  NORTH  AMERICAN  EARTHQUAKES 

Comparison  With  Other  Scaling  Relations 

The  data  used  in  this  study  suggest  several  differences  to  the  revised  mid-plate  scal¬ 
ing  relation  of  Nuttli  et  al.  (1989).  Fault  lengths  predicted  by  Nuttli  et  al.  (1989)  tend  to 
be  shorter  (Figure  12),  fault  areas  smaller  (Figure  8),  and  average  displacements  larger 
(Figure  9)  than  suggested  from  world-wide  data  of  intra-continental  earthquakes.  On  the 
other  hand,  four  spectral  scaling  relations  (Hasegawa,  1983;  Atkinson,  1984;  Nuttli  et  al. , 
1989;  and  this  study)  predict  very  similar  durations  in  the  moment  range  under  study  (Fig¬ 
ure  11),  and  the  data  of  Table  4  seem  to  be  fit  by  any  of  these  relations  (c.f.,  Joyner, 
1984). 

The  data  can  be  interpreted  using  a  constant  stress  drop  model  as  well  as  a  variable 
stress  drop  scaling.  However,  the  spectral  scaling  relations  of  both  models  are  very  simi¬ 
lar.  A  spectral  scaling  with  slopes  between  3.0  -  3.3  can  explain  the  data.  A  spectral  slope 
of  3.0  was  inferred  by  Somerville  et  al.  (1987)  from  an  analysis  of  a  much  smaller  data 
set.  A  spectral  slope  of  3.6  as  inferred  for  earthquakes  in  eastern  North  America 
(Hasegawa,  1983;  Atkinson,  1984;  Nuttli  et  al. ,  1989)  is  not  supported  by  the  data. 

We  notice  however  that  events  in  Table  4  scatter  about  the  mean  of  a  scaling  rela¬ 
tion.  For  example,  events  #  9,  14,  and  39  have  smaller  fault  dimensions,  larger  average 
displacements,  and  larger  stress  drops;  and  hence  are  better  described  by  the  relation  of 
Nuttli  et  al.  (1989).  Events  like  these  caution  the  application  of  scaling  relations  in  earth¬ 
quake  hazard  studies  since  anomalous  events  with  small  fault  dimensions  could  produce 
destructive  earthquakes. 

The  linear  regression  of  observed  body-wave  magnitudes  versus  seismic  moment  has 
a  smaller  slope  than  that  of  Nuttli  et  al.  (1989),  who  did  not  take  the  saturation  of  the 
magnitude  scale  into  account.  The  relation  of  this  study  will  overpredict  magnitudes  for 
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small  seismic  moments,  and  a  straight  line  is  not  suitable  for  fitting  the  data  for  all  seismic 
moments. 

Consistency  Check  on  Scaling  Laws 

Next,  we  will  predict  earthquake  magnitudes  for  events  in  eastern  North  America 
using  the  proposed  scaling  relations  and  compare  these  synthetic  magnitudes  to  the  model 
input  (to  check  for  internal  consistency)  and  to  observations  in  eastern  North  America. 

First,  we  modified  random  process  theory  (e.g.,  (19)  in  Herrmann,  1987)  to  include 
our  proposed  scaling  models.  Following  Herrmann  (1987),  we  chose  /max  =  50  Hz  and 
Q(f)  =  900  *  /a2°.  The  ground  motions  generated  by  random  process  theory  were  used  to 
determine  mLg  following  Herrmann  and  Kijko  (1983).  We  calculated  synthetic  time  his¬ 
tories  at  only  one  epicentral  distance  (300  km)  since  values  are  essentially  indepen¬ 
dent  of  distance  (Boore  and  Atkinson,  1987).  Lg-wave  magnitudes  from  random  process 
theory  using  our  constant  stress  drop  scaling  are  displayed  as  in  Table  5  and  6  for 
constant  and  increasing  stress  drop,  respectively.  The  latter  are  essentially  identical  to 
values  given  by  Toro  and  McGuire  (1987)  for  a  constant  stress  drop  of  100  bar.  This 
equivalence  also  holds  with  respect  to  out  constant  stress  drop  scaling,  except  for  large 
seismic  moments  where  our  magnitudes  are  somewhat  lower.  Comparing  the  synthetic 
Lg-wave  magnitudes  to  the  observed  teleseismic  body-wave  magnitudes,  mb,  shows  that 
the  postulated  equivalence  of  both  magnitude  scales  does  not  hold  for  large  events.  The 
obvious  reason  is  the  saturation  of  the  mb  scale. 

As  a  second  approach,  we  employed  a  deterministic  modeling  technique  for  a  finite 
source  (Herrmann  and  Jost,  1988;  lost  and  Herrmann,  1988;  Jost,  1989).  This  technique 
was  developed  to  study  low-frequency  ground  motions  at  regional  distances  due  to  large 
mid-plate  earthquakes.  We  assumed  that  any  earthquake  is  a  superposition  of  subevents, 
and  that  the  rupture  initiates  from  a  point  that  is  equidistant  from  both  fault  ends  The  rup¬ 
ture  initiation  depth  was  held  constant  at  10  km  as  the  fault  dimensions  change.  The  rup¬ 
ture  initiated  at  the  bottom  of  the  fault  plane  for  small  events.  For  larger  events,  the  point 
of  rupture  initiation  moved  more  to  the  middle  of  the  fault  plane  (the  upper  fault  boundary 
was  fixed  at  3  km  below  earth-surface).  Two  focal  mechanisms  were  considered:  a 
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steeply  dipping  (80°)  strike-slip  and  a  45°  dip-slip  faults.  The  sub-events  on  the  fault 
plane  were  generated  using  normal  mode  theory  (Herrmann  and  Wang,  1985)  and  the 
Central  United  States  earth-model  (Herrmann,  1986).  Assuming  Qa  =  2  Qp,  we  used  Qa 
values  of  100  for  the  first  0.5  km,  500  for  the  second  0.5  km,  1000  for  the  next  19  km, 
and  4000  thereafter.  Greens-functions,  serving  as  subevents,  were  calculated  at  source 
depths  of  4,  6,  8,  10,  20,  30,  and  40  km.  These  source-depth  values  were  associated  with 
the  depths  of  the  subevents  on  the  finite  fault  Rupture  velocities  (90  %  the  shear  wave 
velocity)  and  the  seismic  moments  of  the  subevents  were  somewhat  randomized  to  prevent 
artificial  periodicities  due  to  the  uniform  grid.  Synthetic  time  histories  were  generated  at 
16  azimuths,  equally  distributed  about  the  source.  Lg-wave  magnitudes  (Herrmann  and 
Kijko,  1983)  were  calculated  for  each  trace  and  averaged  with  respect  to  azimuth  and 
focal  mechanism.  Furthermore,  teleseismic  magnitudes  were  determined  using  the  same 
technique  employing  Green’s  functions  at  50°  (Hudson,'  1969).  In  addition,  a  maximum 
mb,  i.e.,  rhb,  was  determined  following  Houston  and  Kanamori  (1986).  The  results  of  this 
magnitude  estimation  are  displayed  in  Tables  5  and  6,  labeled  mffi,  mb3\  and  mb . 

Lg-wave  magnitudes  from  random  process  theory  and  those  from  the  finite  source 
modeling  essentially  agree.  Observed  and  predicted  teleseismic  body-wave  magnitudes 
also  show  reasonable  agreement  thus  supporting  the  internal  consistency  of  the  scaling. 
Values  of  mb  and  mLg  are  not  similar.  However  values  of  rhb  and  mLg  are  comparable  in 
the  full  moment  range  under  study. 

The  results  of  our  kinematic  source  modeling  showed  that  the  Lg-wave  magnitudes 
of  a  strike-slip  fault  are  0.3  magnitude  units  smaller  than  those  of  a  45°  dip-slip  fault.  At 
teleseismic  distances,  the  difference  in  magnitudes  between  strike-slip  and  dip-slip  faulting 
was  found  to  be  somewhat  larger,  the  strike-slip  mb  being  smaller  by  0.5  magnitude  units. 
The  standard  deviations  in  magnitude  values  in  Tables  5  and  6  reflect  the  focal  mechanism 
and  azimuthal  differences  in  the  observations.  They  are  somewhat  large  since  our  numeri¬ 
cal  technique  does  not  consider  scattering  which  may  smooth  the  radiation  pattern.  Hence 
we  observed  very  pronounced  nodes,  which  arc  not  present  in  observations,  at  certain 
azimuths  for  the  strike  slip  case. 
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Comparison  to  Eastern  North  American  Earthquakes 

Next,  we  will  compare  the  magnitude  relation  in  Table  5  to  well  observed  events  in 
eastern  North  America.  The  Oct.  21,  1965  Missouri  event  had  mb  =  4.85  ±  0.23,  = 

5.04  ±  0.12,  Ms  =  4.13  ±  0.32,  and  M0  =  9.0  *  10 22  dyne-cm  (Nuttli,  1973b,  1983). 
Accepting  the  seismic  moment  as  error  free,  we  note  that  our  relations  predict  somewhat 
larger  magnitudes  than  observed  for  this  event,  where  predictions  in  Table  6  fit  the  data 
slightly  better  than  those  in  Table  5.  The  Nov.  9,  1968  Illinois  event  had  mb  =  5.50  ± 
0.40,  mLg  =  5.38  ±  0.23.  Ms  =  5.26  ±  0.28,  and  Af  0  =  9.7  *  1023  dyne-cm  (Nuttli,  1973b, 
1983).  These  observed  source  parameters  are  in  agreement  with  those  in  Tables  5  and  6, 
except  the  Lg-wave  magnitude,  which  is  overpredicted  by  0.5  magnitude  units  in  both 
tables.  The  Jan.  9,  1982  New  Brunswick  event  had  mb  =  5.8,  mLg  =  5.7,  Ms  =  5.1,  and 
Mq  =  2.2  *  1 024  dyne-cm  (North  et  al. ,  written  communication,  1989).  For  this  event,  all 
source  parameters  essentially  agree  with  the  predictions  of  Tables  5  and  6.  Finally,  the 
Nov.  25,  1988  Saguenay  event  had  mb  =  5.9,  mLg  =  6.5,  Ms  =  5.7,  and  M0  =  8.0  *  1024 
dyne-cm  (North  et  al.,  written  communication,  1989).  We  note  a  reasonable  agreement 
with  Tables  5  and  6;  however,  the  observed  Lg-wave  magnitude  of  this  event  is  larger 
than  predicted  by  the  scaling  relations  (North  et  al.,  written  communication,  1989).  In 
conclusion,  the  scaling  relations  proposed  seem  to  predict  magnitudes  for  eastern  North 
America  reasonably  well,  with  a  tendency  though  to  overestimate  mLg  by  0.3  magnitude 
units  for  the  central  United  States  and  to  underestimate  it  for  Saguenay. 

Using  the  proposed  scaling  relations  for  intr>-continental  earthquakes  (Table  5),  we 
can  estimate  source  parameters  for  some  of  the  largest  historical  events  in  the  eastern 
United  States.  However,  published  seismic  moments  for  these  events  cannot  be  used  since 
they  have  been  estimated  from  magnitudes  using  empirical  relations.  The  observations 
available  for  these  historical  events  are  intensity  values  which  have  been  used  to  deter¬ 
mine  magnitudes  such  as  mb,  mLg ,  or  Ms  (e.g.,  Nuttli,  1973a;  Nuttli  et  al. ,  1979). 

For  the  1886  earthquake  near  Charleston,  S.  C.,  Nuttli  (1983)  estimated  a  Lg-wave 
magnitude  of  6.7  (Campbell,  1986).  This  magnitude  is  a  slight  revision  of  his  earlier  esti¬ 
mate  of  6.6  based  on  a  relation  of  intensity  versus  magnitude  (Nuttli  et  al. ,  1979).  Using 
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Table  5,  we  would  expect  a  body-wave  magnitude  of  6.0,  an  Ms  of  6.8,  a  seismic 
moment  of  1.0  *  1026  dyne-cm,  a  fault  length  of  23  km,  a  fault  width  of  11  km,  and  an 
average  displacement  on  the  fault  plane  of  1.1  m.  Assuming  that  present  day  micro¬ 
earthquake  activity  occurs  on  the  fault  plane  that  ruptured  in  1886,  we  can  compare  our 
prediction  of  fault  dimensions  with  observations  by  Shedlock  (1987).  She  inferred  a  length 
extent  of  22  km  and  depth  extent  of  12  km  for  the  hypocenter  distribution  of  recent 
micro-earthquakes  in  the  hypocenter  region  of  the  1886  event  (c.f.,  Talwani,  1982).  Ear¬ 
lier  estimates  by  Nuttli  et  al.  (1979)  and  Nuttli  (1983)  indicated  a  fault  length  of  30  km 
(Nuttli  et  al .,  1989).  We  conclude  that  fault  dimensions  from  recent  micro-earthquakes 
are  in  good  agreement  to  predictions  from  our  proposed  scaling  model. 

Next,  we  address  the  earthquakes  near  New  Madrid  in  1811-1812.  Based  on  a  rela¬ 
tion  between  intensity  and  Lg-wave  magnitude,  Nuttli  et  al.  (1979)  estimated  a  Lg-wave 
magnitude  of  7.3  for  the  earthquake  on  December  16,  1811,  7.2  for  the  event  on  January 
23,  1812,  and  7.5  for  the  earthquake  on  February  7,  1812.  These  estimates  were  based  on 
a  curve  of  intensity  falloff  versus  distance  calibrated  by  using  intra-  and  interplate  earth¬ 
quakes  (Nuttli  et  al. ,  1979).  Assuming  that  including  events  from  California  for  the  cali¬ 
bration  procedure  in  Nuttli  et  al.  (1979)  did  not  bias  the  magnitude  estimates  for  New 
Madrid,  we  estimate  the  following  source  parameters  for  the  historical  events  (Table  5). 
For  the  Dec.  16  earthquake,  we  obtain  an  mb  of  6.2,  an  Ms  of  7.6,  a  seismic  moment  1.0 
*  1027  dyne-cm,  a  fault  length  of  62  km,  a  fault  width  of  22  km,  and  a  mean  displacement 
on  the  fault  plane  of  2.3  m.  A  similar  estimate  would  also  hold  if  we  assume  increasing 
stress  drop  scaling  (Table  6).  Noting  the  slight  differences  between  Lg-wave  magnitudes 
estimated  from  finite  source  modeling  and  random  process  theory,  we  get  somewhat  larger 
source  parameters  by  assuming  the  magnitude  from  tandom  process  theory.  The  January 
23,  1812  earthquake  (Nuttli  et  al.,  1979)  was  only  slightly  smaller  than  the  shock  in 
1811,  and  had  probably  very  similar  source  parameters.  However,  the  event  on  Feb.  7, 
1812  was  larger.  From  Table  5,  we  get  mb  =  6.3,  Ms  =  8.0,  M0  =  3.0  *  1027  dyne-cm, 
fault  length  =  100  km,  fault  width  =  32  km,  and  mean  displacement  on  the  fault  =  3.4  m. 
These  estimates  could  also  be  somewhat  larger  if  Lg-wave  magnitudes  from  random  pro¬ 
cess  theory  rather  than  from  finite  source  modeling  would  be  used. 
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Studies  of  microseismicity  in  the  New  Madrid  Seismic  Zone  (Stauder  et  al. ,  1976) 
indicate  a  length  of  about  125  km  for  the  southern  segment  which  is  associated  with  the 
event  on  December  16,  1811  and  its  two  major  aftershocks.  This  segment  of  the  fault 
could  easily  contain  a  rupture  length  of  60  -100  km  corresponding  to  an  event  with 
seismic  moment  between  1.0  *  1027  dyne-cm  and  3.0  *  1027  dyne-cm.  The  whole  length  of 
125  km  could  have  ruptured  in  several  events  rather  than  in  one.  The  above  estimate  of 
fault  width  is  more  problematic  since  present  day  seismicity  in  the  New  Madrid  Seismic 
Zone  extents  only  to  a  depth  of  about  14  km  (Himes  et  al. ,  1988).  On  the  other  hand, 
the  November  9,  1968  Illinois  earthquake  had  a  confirmed  depth  of  25  km,  suggesting  that 
deeper  parts  in  the  mid-continental  crust  can  be  seismogenic.  The  earthquake  near  New 
Madrid  on  January  23,  1812  is  associated  with  the  central  section  of  the  New  Madrid 
Seismic  Zone.  Its  length  is  about  60  km  which  agrees  well  with  our  above  estimate.  For 
the  third  event  on  February  7,  1812,  the  length  of  microseismicity  of  about  85  km  on  the 
northern  part  of  the  New  Madrid  Seismic  Zone  (Nuttli,  1983)  is  only  slightly  smaller  than 
our  estimate  from  above.  Considering  the  variability  in  our  source  parameters,  even  a 
seismic  moment  of  of  8  x  1027  dyne-cm  as  proposed  by  Johnston  (1989)  and  Nuttli  et  al. 
(1989)  could  be  reasonable  for  this  event. 

These  estimates  of  source  parameters  apply  only  if  the  corresponding  events  were 
average  intra-continental  earthquakes.  As  discussed  before,  actual  events  show  a  scatter  of 
source  parameter  values  such  that  each  estimate  of  source  dimension,  average  displace¬ 
ment,  and  stress  drop  could  vary  by  a  factor  of  2  (Bonilla  et  al. ,  1984).  Another  problem 
is  related  to  the  thickness  of  the  seismogenic  zone  which  can  reach  up  to  45  km  according 
to  the  intra-continental  scaling  relations  assuming  constant  stress  drop  (Table  5).  In 
regions  like  the  central  United  States,  this  value  for  the  thickness  of  the  seismogenic  zone 
seems  too  large  (Joyner,  1984).  It  is  interesting  to  note  that  the  increasing  stress  drop 
scaling  for  intra-continental  earthquakes  (Table  6)  would  suggest  fault  widths  of  10,  17, 
and  23  km  for  the  events  in  1886,  1811,  and  1812  (February  7),  respectively. 

Purcaru  and  Berckhemer  (1982)  expressed  doubts  that  a  classification  of  earthquakes 
according  to  plate  tectonics  (e.g.,  mid-plates,  plate  margins)  leads  to  meaningful  results. 
This  contrasts  with  other  studies  suggesting  a  grouping  according  to  type  of  plate 


boundary  (e.g.,  Kanamori  and  Anderson,  1975;  Singh  et  al.,  1980;  Nuttli,  1983,  1985). 
As  more  and  better  quality  data  become  available,  scaling  relations  will  need  further 
updating  to  investigate  possible  regional  differences  in  scaling  (Acharya,  1979;  Hasegawa, 
1983;  Bonilla  et  al. ,  1984).  Further  refinement  could  also  come  from  including  focal 
mechanisms  or  other  parameters  into  scaling  relations  (Purcaru  and  Berckhemer,  1982; 
Bonilla  et  al. ,  1984). 


CONCLUSION 

A  data  set  of  59  intra-continental  earthquakes  has  been  used  to  derive  seismic  scaling  rela¬ 
tions.  The  data  can  be  explained  by  a  constant  stress  drop  scaling  (slope  of  3.0)  as  sug¬ 
gested  by  Somerville  et  al.  (1987).  On  the  other  hand,  an  increasing  stress  drop  scaling 
(slope  of  3.3)  cannot  be  ruled  out.  This  range  of  variability  in  spectral  scaling  has  also 
been  observed  for  the  Nahanni  earthquake  sequence  (Boore  and  Atkinson,  1989)  assuming 
these  events  to  be  intra-continental.  Comer  periods  are  very  similar  in  this  range  of  spec¬ 
tral  scaling  exponents.  However,  the  world-wide  data  do  not  support  a  spectral  slope  of 
3.6  as  suggested  for  eastern  North  America  by  Nuttli  et  al.  (1989).  Comparing  Table  5 
to  Table  1  in  Nuttli  (1985)  for  plate  margin  events  indicates  smaller  fault  lengths  for 
intra-continental  versus  plate-margin  events.  On  the  other  hand,  stress  drops  and  average 
displacements  do  not  seem  to  be  fundamentally  different. 

The  proposed  scaling  relations  have  been  used  to  determine  synthetic  Lg-wave  and 
teleseismic  body-wave  magnitudes  to  check  for  internal  consistency.  Observed  magni¬ 
tudes  of  events  from  North  America  show  agreement  with  predictions  of  the  scaling  rela¬ 
tions  considered  the  scatter  in  source  data.  Source  parameters  also  have  been  derived  for 
the  largest  historical  earthquakes  in  the  central  and  eastern  United  States  using  the  pro¬ 
posed  scaling  relations.  Fault  rupture  lengths  show  a  reasonable  agreement  with  observed 
distributions  of  microseismicity.  The  relations  based  on  increasing  stress  drop  predict 
fault  rupture  widths  that  are  in  closer  agreement  with  observed  distributions  of 
microseismicity  in  the  central  and  eastern  United  States  than  those  of  the  constant  stress 
drop  model. 
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AVERAGE  SOURCE  PARAMETERS  FOR  MID-PLATE 
EARTHQUAKES  (NUTTLI  ET  AL.,  1989)* 


5 


average  fault  displacement  on  the  rupture  surface. 


TABLE  2 


SIMILARITY  CONDITIONS  FOR  MID-PLATE  EARTHQUAKES 


log  Mo 

[dyne-cm] 

tr** 

[s] 

W/L** 

u/L** 

io-5 

vrIj/L** 

1.00 

1.4 

0.34 

1.00 

2.4 

0.38 

1.00 

4.3 

0.35 

0.19 

1.00 

7.5 

0.37 

0.37 

0.97 

9.7 

0.36 

0.71 

0.73 

11.5 

0.32 

1.38 

0.54 

13.8 

0.27 

26.5 

1.87 

0.48 

16.5 

0.26 

2.56 

0.41 

18.8 

0.24 

27.5 

3.55 

0.35 

19.6 

0.22 

28.0 

5.01 

0.32 

22.4 

0.21 

**The  meaning  of  the  symbols  is: 
tj  Rise  Time  (Circular  Fault,  Geller,  1976) 

W/L  Aspect  Ratio  (Kanamori  and  Anderson,  1975) 
u/L  Strain  Drop  (Kanamori  and  Anderson,  1975) 

vR  t/L*  Dynamic  Similarity  (Kanamori  and  Anderson,  1975) 

*  Ac^umes  a  rupture  velocity  of  3.15  km/s. 


TABLE  3 


HYPOCENTERS  AND  FOCAL  MECHANISMS 
FOR  INTRA-CONTINENTAL  EARTHQUAKES 


No. 

Region 

1 

Bolnai 

2 

Bolnai 

3 

Kebin 

4 

Haiyuan 

5 

Ku-long 

6 

Fuh-Yun 

7 

Khait 

8 

S.  Tibet 

9 

Mato  Grosso 

10 

Muya 

11 

Gobi- Altai 

12 

Mato  Grosso 

13 

Hsingtai 

14 

Hsingtai 

15 

Hsingtai 

16 

Hsingtai 

17 

Mogod 

18 

Garni 

19 

Ganzi 

20 

Koyna 

21 

Acre 

22 

Meckering 

23 

Tien  Shan 

24 

Bohai 

25 

Ceres 

26 

Tonghai 

27 

Lake  McKay 

28 

Alma-ata 

29 

Tien  Shan 

30 

Dzhanbul 

31 

Tien  Shan 

32 

Simpson  Desert 

33 

Luhuo 

34 

Luhuo 

35 

Tien  Shan 

36 

Yunnan 

37 

Markansu 

38 

Colombia 

39 

Haicheng 

40 

Gazli 

41 

Gazli 

42 

Lung  Ling 

43 

Lung  Ling 

44 

Tangshan 

45 

Luanxian 

Date  Time 
D/M/Y 


Lat. 


Long.  h  T  -  axis  P  -  axis 
[km]  TRD  PLG  TRD  PLG 


09/07/05 

23/07/05 

03/01/11 

16/12/20 

22/05/27 

10/08/31 

10/07/49 

18/11/51 

31/01/55 

27/06/57 

04/12/57 

13/02/64 

07/03/66 

22/03/66 

22/03/66 

26/03/66 

05/01/67 

30/08/67 

30/08/67 

10/12/67 

27/08/68 

14/10/68 

11/02/69 

18/07/69 

29/09/69 

04/01/70 

24/03/70 

05/06/70 

23/03/71 

10/05/71 

09/04/72 

28/08/72 

06/02/73 

07/02/73 

02/06/73 

10/05/74 

11/08/74 

27/09/74 

04/02/75 

08/04/76 

17/05/76 

29/05/76 

29/05/76 

27/07/76 

28/07/76 


09:40:24.0 

02:46:24.0 

23:25:45.0 

12:05:48.0 

22:32:42.0 

21:18:40.0 

03:53:36.0 

09:35:47.0 

05:03:07.0 

00:09:28.0 

03:37:45.0 

11:21:46.0 

21:29:14.0 

08:11:36.0 

08:19:46.0 

15:19:04.0 

00:14:40.1 

04:22:05.1 

11:08:50.1 

22:51:24.3 

05:17:36.0 

02:58:51.8 

22:08:51.0 

05:24:45.0 

20:02:32.1 

17:00:39.0 

10:35:16.9 

04:53:07.4 

20:47:16.0 

14:51:45.0 

04:10:48.9 

02:18:59.4 

10:37:10.1 

16:06:25.0 

23:57:02.4 

19:25:17.0 

01:13:55.0 

04:09:02.0 

11:36:07.5 

02:40:25.0 

02:58:40.0 

12:23:18.4 

14:00:19.4 

19:42:54.6 

10:45:37.2 


49.00N 
49.00N 
42.80N 
36.62N 
38.05N 
46.89N 
39.27N 
30.98N 
12.42S 
56.20N 
45.31N 
18.06S 
37.35N 
37.50N 
37.53N 
37.68N 
48.15N 
31.60N 
31.70N 
17.38N 
8.90S 
31.54S 
41.42N 
38.43N 
33.16S 
24. 10N 
22.08S 
42.48N 
41.42N 
42.85N 
42.09N 
25.01S 
31.40N 
31.50N 
44.14N 
28.19N 
39.38N 
2.72N 
40.65N 
40.36N 
40.37N 
24.5  IN 
24.54N 
39.60N 
39.7  IN 


98.00E 

19 

143 

94.50E 

19 

143 

77.30E 

25a 

0 

105.40E 

25a 

306 

102.37E 

25a 

316 

90.06E 

25a 

115 

70.59E 

25a 

0 

91.49E 

25a 

233 

57.30W 

33h 

342 

116.59E 

25a 

173 

99.21E 

25a 

312 

56.69W 

5h 

345 

114.92E 

10G 

343 

115.08E 

9g 

163 

115.05E 

9g 

340 

115.27E 

15g 

161 

102.90E 

241 

135 

100.30E 

8y 

141 

100.30E 

10Y 

323 

73.75E 

4k 

68 

72.89W 

26h 

212 

117.00E 

3l 

121 

79.24E 

10" 

205 

119.47E 

10c 

19.31E 

11K 

102.50E 

15b 

71 

126.65E 

8l 

345 

78.71E 

17* 

242 

79.20E 

lla 

44 

71.29E 

15a 

237 

84.58E 

13* 

12 

136.37E 

8l 

325 

100.60E 

10Y 

170 

100.30E 

10Y 

300 

83.60E 

26* 

162 

103.98E 

17c 

73.81E 

23x 

84 

71.37W 

6h 

195 

122.80E 

12n 

336 

63.73E 

10° 

265 

63.44E 

13v 

257 

98.95E 

31 

66 

98.60E 

41 

215 

118.00E 

16v 

170 

118.37E 

16v 

356 

14 

237 

14> 

14 

237 

14> 

90 

160 

0A 

59 

50 

8a 

59 

60 

8a 

0 

25 

0A 

90 

360 

0A 

8 

129 

59a 

86 

142 

4h 

1 

82 

43a 

46 

50 

7f 

46 

254 

1“ 

1 

73 

2° 

4 

73 

7° 

4 

249 

20° 

4 

253 

27° 

0 

45 

0F 

2 

238 

76y 

5 

143 

85y 

2 

336 

36J 

43 

107 

15” 

71 

273 

I7l 

84 

339 

4a 

3 

341 

3b 

83 

78 

1L 

74 

354 

6a 

87 

160 

la 

66 

330 

2a 

86 

190 

4a 

75 

145 

15l 

3 

260 

-Y 

1 

15 

120 

75  Y 

60 

27 

22a 

0 

174 

0X 

2 

285 

5" 

4 

244 

21n 

82 

14 

3t 

75 

124 

10T 

7 

335 

13d 

51 

321 

13d 

16 

260 

2Z 

5 

120 

82z 

77 


46 

Isfara 

31/01/77 

14:26:15.1 

40.1  IN 

70.86E 

12* 

116 

82 

344 

6* 

47 

Gazli 

04/06/78 

19:30:20.4 

40.4  IN 

63.60E 

9q 

103 

81 

343 

5° 

48 

Tien  Shan 

29/03/79 

02:01:32.1 

41.95N 

83.38E 

13* 

318 

82 

164 

8“ 

49 

Cadoux 

02/06/79 

09:47:58.7 

30.73S 

117.21E 

3l 

233 

78 

75 

11L 

50 

W.  Amazon 

06/03/80 

09:46:18.0 

6.17S 

71.16W 

18h 

182 

85 

40 

4h 

51 

El  Asnam 

10/10/80 

12:25:25.1 

36.17N 

1.41E 

14f 

107 

80 

320 

98 

52 

S.  Para 

12/11/80 

21:23:05.0 

8.07S 

50.24W 

33h 

197 

10 

93 

55  H 

53 

Ceara 

20/11/80 

03:29:42.0 

4.30S 

38.40W 

5h 

199 

0 

109 

3h 

54 

Daofu 

23/01/81 

21:13:51.7 

30.93N 

101. 10E 

8y 

6 

0 

96 

0Y 

55 

Paraguay 

08/04/82 

05:58:52.0 

24. 80S 

58. 10W 

12h 

351 

1 

81 

3h 

56 

Sokh 

06/05/82 

15:42:22.2 

40.15N 

71.54E 

20* 

139 

80 

337 

9* 

57 

Codajas 

05/08/83 

06:21:42.0 

3.59S 

62.17W 

23h 

264 

62 

14 

10H 

58 

Guinea 

22/12/83 

04:11:29.2 

11.87N 

13.53W 

ir 

34 

2 

125 

351 

59 

Gazli 

19/03/84 

20:28:38.3 

40.30N 

63.30E 

9q 

346 

83 

125 

5<> 

60 

Marryat  Creek 

30/03/86 

08:53:52.0 

26.30S 

132.77E 

2l 

275 

78 

66 

10L 

61 

Tennant  Creek 

22/01/88 

00:35:58.0 

19.85S 

133.80E 

5P 

177 

73 

20 

16’ 

62 

Tennant  Creek 

22/01/88 

03:57:25.2 

19.80S 

133.91E 

5P 

199 

67 

360 

22' 

63 

Tennant  Creek 

22/01/88 

12:04:57.8 

19.83S 

133.88E 

5P 

177 

73 

20 

16‘ 

References  :  A  Chen  and  Molnar  (1977);  B  Singh  and  Havskov  (1980);  C  Deng  et  al.  (1984);  D  Deng 
et  al.  (1986);  E  Zhang  et  al.  (1987);  F  Okal  (1976);  G  Chung  and  Cipar  (1983);  H  Assumpcao  and 
Suarez  (1988);  I  ISC;  J  Langston  (1976);  K  Somerville  et  al.  (1987);  L  Fredrich  et  al.  (1988);  M 
Vogfjord  and  Langston  (1987);  N  Cipar  (1979);  P  PDE;  Q  Eyidogan  et  al.  (1985);  R  Hartzell  (1980);  T 
Kristy  et  al.  (1980);  U  Liu  and  Kanamori  (1980);  V  Butler  et  al.  (1979);  W  Kanamori  and  Allen  (1986); 
X  Jackson  et  al.  (1979);  Y  Zhou  et  al.  (1983b);  Z  Nabelek  et  al.  (1987);  a  Nelson  et  al.  (1987);  b 
Zhou  et  al.  (1983a);  c  Purcaru  and  Berckhemer  (1982);  d  Okal  and  Stewart  (1981);  e  Langer  et  al. 
(1987);  f  Ck'temas  et  al.  (1982);  g  Deschamps  et  al.  (1982);  h  Ouyed  et  al.  (1981);  i  Chung  et  al. 
(1988);  j  Okal  (1977). 


TABLE  4 


SOURCE  PARAMETERS  FOR  INTRA-CONTINENTAL  EARTHQUAKES 


No.  Date 
D/M/Y 


1 

09/07/05 

2 

23/07/05 

3 

03/01/11 

4 

16/12/20 

5 

22/05/27 

6 

10/08/31 

7 

10/07/49 

8 

18/11/51 

9 

31/01/55 

10 

27/06/57 

11 

04/12/57 

12 

13/02/64 

13 

07/03/66 

14 

22/03/66 

15 

22/03/66 

16 

26/03/66 

17 

05/01/67 

18 

30/08/67 

19 

30/08/67 

20 

10/12/67 

21 

27/08/68 

22 

14/10/68 

23 

11/02/69 

24 

18/07/69 

25 

29/09/69 

26 

04/01/70 

27 

24/03/70 

28 

05/06/70 

29 

23/03/71 

30 

10/05/71 

31 

09/04/72 

32 

28/08/72 

33 

06/02/73 

34 

07/02/73 

35 

02/06/73 

36 

10/05/74 

37 

11/08/74 

38 

27/09/74 

39 

04/02/75 

40 

08/04/76 

41 

17/05/76 

42 

29/05/76 

43 

29/05/76 

44 

27/07/76 

45 

28/07/76 

46 

31/01/77 

mb  Ms 


7.90J 

8.25J 

8.40® 

8.50® 

7.90® 

7.90® 

7.60® 


6.20h 

5.50h 

7.90® 

8.10® 

5.40h 

4.50h 

5.60° 

6.80g 

5.60g 

6.70° 

5.90g 

7.20° 

5.20g 

6.20g 

6.101 

7.40® 

6.101 

6.10y 

5.201 

5.10y 

5.901 

6.501 

4.90h 

3.90h 

5.901 

6.801 

5.801 

6.601 

6.201 

7.30c 

5.601 

6.301 

5.90c 

7.50b 

6.201 

5.901 

5.901 

6.601 

5.801 

5.801 

5.601 

5.401 

5.801 

5.601 

5.301 

6.10p 

7.50y 

5.80p 

5.90y 

5.701 

5.601 

5.801 

6.801 

6.40x 

7.30x 

5.50H 

5.80® 

6.40* 

7.40n 

6.20® 

7.00® 

6.20® 

7.00® 

5.901 

6.901 

5.701 

7.001 

6.30v 

7.70v 

6.101 

7.20v 

6.001 

5.901 

log  Mq 
[dyne-cm] 

T 

[sec] 

L 

[km] 

W* 

[km] 

u 

[m] 

Act 

[bar] 

28.74) 

200) 

50c 

14.0° 

65c 

28.70s 

300) 

50c 

8.2C 

37c 

27.69A 

180A 

40a 

2.3a 

28.08c 

220° 

20° 

8.0e 

55c 

27.63a 

150a 

40a 

2.4a 

31c 

27.93a 

300a 

50a 

1.9a 

30c 

27.38a 

70  A 

30A 

3.7a 

33" 

27.66a 

200a 

27.15a 

35a 

30A 

4.5a 

76c 

28.1  lA 

27 0A 

50A 

3.2a 

26F 

26.00° 

5.0° 

14° 

1.6° 

112° 

25.48° 

4.5° 

13° 

0.7° 

53° 

26.25° 

5.0° 

14° 

2.9° 

194° 

25.15° 

3.5° 

10° 

0.6° 

53° 

26.58a 

40a 

30a 

1.0A 

24c 

25.65y 

2.5Y 

9y 

1.5Y 

35y 

24.34y 

2.5y 

25.5 1J 

6.4® 

32c 

llc 

0.6c 

20k 

26.02l 

5.4l 

37l 

IqM 

2.0l 

100M 

25.30® 

6.9s 

26.90c 

45c 

25c 

2.2c 

25.65k 

5.0k 

70K 

26.94b 

50.0b 

90b 

15b 

2.1b 

27b 

25.07L 

2.5® 

6U 

175u 

25.49a 

3.8s 

24.77a 

7.0s 

24.11s 

1.3s 

24.11s 

1.0s 

24.00L 

1.8l 

27.26y 

32.0y 

90Y 

15y 

3.8y 

50y 

24.77y 

2.5y 

24.32s 

0.8s 

25.8  lc 

45c 

20c 

0.2C 

6C 

26.70x 

30x 

20x 

2.5X 

21x 

26.48N 

7.0n 

22n 

12n 

2.8N 

53n 

26.32*2 

5.0k 

30c 

15c 

400k 

26.30*2 

7.8® 

15® 

10® 

3.3® 

200® 

25.70c 

35c 

20c 

25.98c 

45c 

20c 

27.26v 

25. 02 

80w 

15w 

100® 

26.90v 

10.02 

35v 

15v 

40c 

24.72s 

1.0s 
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47 

04/06/78 

5.90q 

5.70q 

0.2^ 

48 

29/03/79 

5.801 

5.801 

24.36s 

1.3s 

49 

02/06/79 

5.901 

6.001 

25.17L 

6.4l 

15l 

6l 

1.5l 

200k 

50 

06/03/80 

4.80” 

3.60“ 

51 

10/10/80 

6.50f 

7.30f 

26.70s 

12.011 

36s 

10s 

4.0s 

110s 

52 

12/11/80 

4.70“ 

3.30h 

53 

20/11/80 

5.20h 

4.40“ 

23.85“ 

1.1“ 

3“ 

90“ 

54 

23/01/81 

5.701 

6.80y 

26.1 1Y 

12.0y 

44y 

10Y 

1.0Y 

20y 

55 

08/04/82 

4.90h 

4.00h 

56 

06/05/82 

5.501 

5.701 

24.30s 

3.6s 

57 

05/08/83 

5.50“ 

4.50“ 

24.00“ 

1.0“ 

3“ 

100“ 

58 

22/12/83 

6.40e 

6.20e 

25.32k 

5.0* 

9e 

0.1e 

30k 

59 

19/03/84 

6.50q 

7.00q 

26.40^ 

12.0q 

30k 

60 

30/03/86 

5.701 

5.801 

24.76l 

4# 

13l 

3l 

o 

w 

r 

61 

22/01/88 

frlO1* 

6.30p 

25.211 

12* 

62 

22/01/88 

6.10p 

6.40p 

25.431 

9* 

63 

22/01/88 

6.50p 

6.7(f 

25.90i 

141 

*  estimated.  References  :  A  Chen  and  Molnar  (1977);  B  Singh  and  Havskov  (1980);  C  Deng 
et  al.  (1984);  D  Deng  et  al.  (1986);  E  Zhang  et  al.  (1987);  F  Okal  (1976);  G  Chung  and  Cipar 
(1983);  H  Assumpcao  and  Suarez  (1988);  I  ISC;  J  Langston  (1976);  K  Somerville  et  al.  (1987); 
L  Fredrich  et  al.  (1988);  M  Vogfjord  and  Langston  (1987);  N  Cipar  (1979);  P  PDE;  Q  Eyido- 
gan  et  al.  (1985);  R  Hartzell  (1980);  T  Kristy  et  al.  (1980);  U  Liu  and  Kanamori  (1980);  V 
Butler  et  al.  (1979);  W  Kanamori  and  Allen  (1986);  X  Jackson  et  al.  (1979);  Y  Zhou  et  al. 
(1983b);  Z  Nabelek  et  al.  (1987);  a  Nelson  et  al.  (1987);  b  Zhou  et  al.  (1983a);  c  Purcaru  and 
Berckhemer  (1982);  d  Okal  and  Stewart  (1981);  e  Langer  et  al.  (1987);  f  Cistemas  et  al. 
(1982);  g  Deschamps  et  al.  (1982);  h  Ouyed  et  al.  (1981);  i  Chung  et  al.  (1988);  j  Okal 
(1977). 
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AVERAGE  SOURCE  PARAMETERS  FOR  INTRA-CONTINENTAL 
EARTHQUAKES  ASSUMING  CONSTANT  STRESS  DROP  SCALING* 
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Lg-wave  magnitude  from  numerical  modeling  using  a  finite  source  (average  of  SS  and  45°  DS  faults;  Herrmann  and  Kijko,  1983) 

body-wave  magnitude  from  numerical  modeling  using  a  finite  source  (average  of  SS  and  45°  DS  faults) 

maximum  mb  (Houston  and  Kanamori,  1986)  from  numerical  modeling  using  a  finite  source  (average  of  SS  and  45°  DS  faults) 


RAGE  SOURCE  PARAMETERS  FOR  INTRA-CONTINENTAL 
IQUAKES  ASSUMING  INCREASING  STRESS  DROP  SCALING* 
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TABLE  7 

SIMILARITY  CONDITIONS  FOR  INTRA-CONTINENTAL  EARTHQUAKES 

a)  Constant  Stress  Drop 


log  Mo 

[dyne-cm] 

tj** 

[s] 

W/L** 

u/L** 

10"5 

vR  t/L** 

23.0 

0.21 

1.00 

8.5 

0.51 

24.0 

0.45 

0.88 

7.0 

0.43 

25.0 

0.97 

0.65 

5.6 

0.35 

26.0 

2.08 

0.48 

4.8 

0.28 

26.5 

3.06 

0.42 

4.2 

0.25 

27.0 

4.49 

0.35 

3.7 

0.23 

27.5 

6.59 

0.32 

3.4 

0.21 

28.0 

9.68 

0.27 

3.0 

0.18 

b)  Increasing  Stress  Drop 


log  M0 

[dyne-cm] 

v** 

[s] 

W/L** 

u/L** 

io-5 

vR  t/L** 

23.0 

0.28 

1.00 

4.0 

0.59 

.  "Ill 

0.55 

1.00 

4.5 

0.53 

1.10 

0.67 

4.3 

0.39 

2.19 

0.43 

4.2 

0.30 

26.5 

3.09 

0.34 

3.9 

0.26 

4.36 

0.27 

3.9 

0.22 

27.5 

6.16 

0.23 

3.8 

0.19 

28.0 

8.70 

0.18 

3.7 

0.17 

**The  meaning  of  the  symbols  is: 
tj  Rise  Time  (Circular  Fault,  Geller,  1976) 

W/L  Aspect  Ratio  (Kanamori  and  Anderson,  1975) 
u/L  Strain  Drop  (Kanamori  and  Anderson,  1975) 

vR  t/L*  Dynamic  Similarity  (Kanamori  and  Anderson,  1975) 
*  Assumes  a  rupture  velocity  of  3.15  Icm/s. 
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Figure  1.  Reduced  far-field  displacement  spectra  for  central  and  eastern  North  American 
earthquakes  defining  the  revised  mid-plate  scaling  relation  by  Nuttli  et  al.  (1989). 
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Figure  2.  Relation  between  the  logarithm  of  seismic  moment,  M0,  and  comer  periods,  Tc , 
of  the  revised  mid-plate  scaling  relation  by  Nuttli  et  al.  (1989).  The  solid  line  is  for  the 
comer  period  T 02  (equations  la,  3a,  6a),  long  dashes  are  for  T01  (equations  la,  2,  5),  and 
short  dashes  are  for  Tl2  (equations  la,  4,  7),  where  the  indices  refer  to  the  slopes  of  the 
reduced  far-field  displacement  spectra  (Figure  1). 
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Figure  3.  Source  duration  versus  seismic  moment  for  eastern  North  American  earth¬ 
quakes.  Data  are  from  Somerville  et  al.  (1987,  circles  for  events  with  stress  drops  larger 
or  equal  200  bar,  x’s  for  events  with  stress  drops  smaller  than  200  bar),  Fletcher  (1982, 
triangles),  and  Shin  and  Herrmann  (1987,  crosses).  Included  are  the  relations  by 
Hasegawa  (1983,  long  dashes),  Nuttli  (1983,  intermediate  dashes),  Somerville  et  al. 
(1987,  solid  line),  and  Nuttli  et  al.  (1989,  short  dashes). 
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Figure  4.  Static  stress  drop  versus  source  duration  (i.e.  comer  period  T ^  for  eastern 
North  American  earthquakes.  Data  are  from  Somerville  et  al.  (1987,  circles  for  events 
with  stress  drops  larger  or  equal  200  bar,  x’s  for  events  with  stress  drops  smaller  than 
200  bar),  and  Shin  and  Herrmann  (1987,  crosses).  Included  are  the  relations  by  Somer¬ 
ville  et  al.  (1987,  solid  line)  and  Nuttli  et  al.  (1989,  dashes). 
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Figure  6.  Observed  P-wave  magnitudes,  mb,  versus  seismic  moment  from  Table  4.  The 
solid  line  is  for  the  regression  of  mb  versus  log  M 0,  (22a),  the  long  dashes  for  the  regres¬ 
sion  of  log  M0  versus  mb,  (22b).  Short  dashes  are  from  the  revised  mid-plate  scaling  rela¬ 
tion  by  Nuttli  et  al.  (1989),  (23). 
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Figure  7.  Observed  surface-wave  magnitudes,  Ms,  versus  seismic  moment  from  Table  4. 
The  solid  line  is  from  the  least-squares  fit,  (25),  short  dashes  are  from  the  revised  mid¬ 
plate  scaling  relation  by  Nuttli  et  al.  (1989),  (13),  and  long  dashes  are  for  Mw  (Singh 
and  Havskov,  1980),  (21). 
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Figure  8.  Fault  areas,  A,  versus  seismic  moment  from  Table  4.  The  solid  line  is  from  the 
least-squares  fit,  (26),  dash-dots  are  from  relations  assuming  constant  stress  drop  scaling, 
(28),  long  dashes  from  relations  assuming  increasing  stress  drop,  (30).  Short  dashes  are 
from  the  revised  mid-plate  scaling  relation  by  Nuttli  et  al.  (1989). 
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Figure  9.  Average  displacements,  u ,  versus  seismic  moment  from  Table  4.  The  solid  line 
is  from  the  least-squares  fit,  (27),  dash-dots  are  from  relations  assuming  constant  stress 
drop  scaling,  (29),  long  dashes  from  relations  assuming  increasing  stress  drop,  (31).  Short 
dashes  are  from  die  revised  mid-plate  scaling  relation  by  Nuttli  et  al.  (1989). 
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Figure  10.  Observed  static  stress  drops,  A o,  versus  seismic  moment  from  Table  4.  For 
intra-continental  earthquakes,  dash-dots  indicate  constant  stress  drop  (45  bar,  (32)),  long 
dashes  increasing  stress  drop,  (33).  Short  dashes  are  from  the  revised  mid-plate  scaling 
relation  by  Nuttli  et  al.  (1989),  and  the  solid  line  is  from  Somerville  et  al.  (1987). 
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Figure  11.  Observed  durations,  t,  versus  seismic  moment  from  Table  4.  Dash-dots  are 
from  relations  assuming  constant  stress  drop  scaling,  (36),  long  dashes  from  relations 
assuming  increasing  stress  drop,  (37).  Short  dashes  are  from  the  revised  mid-plate  scaling 
relation  by  Nuttli  et  al.  (1989),  and  the  solid  line  is  from  Somerville  et  al.  (1987). 


LOG ( L)  [cm] 


n 

9 


2  2  23  24  25  26  27  2  8  29  30 

LOG(M0)  [dyne -cm] 

Figure  12.  Observed  fault  lengths,  L,  versus  seismic  moment  from  Table  4.  The  solid 
line  is  from  the  least-squares  fit,  (40),  and  short  dashes  are  from  the  revised  mid-plate 
scaling  relation  by  Nuttli  et  al.  (1989). 
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Figure  13.  Comer  periods  versus  seismic  moments  derived  from  scaling  relations  for  intra- 
continental  earthquakes,  (a)  Constant  stress  drop  scaling:  The  solid  lines  are  from  top  to  bottom  for 
7oi.  (41),  T0 2,  (36),  and  T 12,  (43),  respectively.  For  comparison,  short  dashes  are  from  the  revised 
mid-plate  scaling  relation  by  Nutdi  et  at.  (1989),  i.e.,  Figure  2,  for  Tou  T0i,  and  Tn,  respectively, 
Observed  comer  periods  are  displayed  for  events  40  and  41  (Upadhyay  and  Duda,  1980). 
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Figure  13.  Comer  periods  versus  seismic  moments  derived  from  scaling  relations  for  intra¬ 
continental  earthquakes,  (b)  Comparison  of  comer  periods  derived  from  constant  stress  drop  (solid 
lines)  and  increasing  stress  drop  scaling  (short  dashed  lines). 
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Ground  roll:  rejection  using  polarization  filters 


Chiou-Fen  Shieh*  and  Robert  B.  Herrmann^ 

ABSTRACT 

Ground  roll  noise  on  land  data  sets  overwhelms  the  desired  reflection 
seismic  signal  unless  special  steps  are  taken  in  data  acquisition  and  process¬ 
ing  to  control  it.  This  is  usually  done  in  the  field  by  the  design  of  group 
arrays  for  data  acquisition.  On  the  other  hand,  if  multicomponent  data  are 
acquired,  it  is  possible  to  remove  the  ground  roll  during  processing  by  using 
polarization  analysis.  Even  though  this  processing  is  computationally  inten¬ 
sive,  the  potential  exists  for  obtaining  results  similar  to  conventional  data 
acquisition  with  the  deployment  of  fewer  sensors  in  the  field  with  minimal 
group  array  effects  as  well  as  deriving  new  information. 

We  describe  a  two-dimensional  polarization  filter  analysis  for  use  with 
vertical  and  inline  sensors.  A  time-domain  spectral  matrix  technique  is 
developed  to  account  for  the  fact  that  the  recorded  seismic  signal  is  the  super¬ 
position  of  multiple  signals  in  the  time  domain,  each  with  different  frequency 
content  and  time-varying  polarization.  This  technique  is  implemented  by 
decomposing  the  signal  into  individual  frequency  components  by  using  nar¬ 
row  bandpass  filters  and  defining  the  polarization  characteristic  using  sliding 
time  windows.  We  show  that  both  incoherent  noise  and  specific  linearly 
polarized  constituents  can  be  successfully  filtered. 

INTRODUCTION 

Ground  roll  can  be  eliminated  from  a  seismic  trace  by  exploiting  its 
specific  characteristics.  For  example,  one  may  focus  on  the  dispersion,  and 
design  numerical  filters  to  remove  the  dispersed  signal  from  the  trace 
(Herrmann  and  Russell,  1990).  An  alternative  approach  is  to  realize  that  the 
recorded  seismic  wave  field  may  be  represented  to  a  first  approximation  by 
the  superposition  of  surface-wave  modes  (Harvey,  1981).  The  surface-wave 
signal  in  an  isotropic  elastic  medium  has  the  feature  that  the  vertical  and 
inline  components  are  90°  out  of  phase  and  that  the  elliptidty  is  frequency 
and  mode  dependent  (Mooney  and  Bolt,  1966;  Harkrider,  1970).  This  suggests 
that  an  analysis  of  the  instantaneous  particle  motion,  or  its  polarization  if  the 
motion  is  elliptical,  may  lead  to  another  method  of  ground  roll  removal. 

Polarization  analysis  has  been  used  to  introduce  extra  parameters  for  the 
description  of  complex  wave  fields.  (Renfc  et  al,  1986).  If  one  can  calculate 
those  parameters  accurately,  identification  and  isolation  of  signal  components 
can  be  achieved.  Methods  have  been  presented  to  perform  this  analysis  in 
either  the  frequency  or  time  domains. 

•Formerly  Department  of  Earth  and  Atmospheric  Sciences,  Saint  Louis  University; 
currently  Taipei  Teachers  College,  Taipei,  Taiwan,  ROC. 

§Department  of  Earth  and  Atmospheric  Sciences,  Saint  Louis  University,  3507  Laclede 
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A  spectral  matrix  method  was  first  introduced  by  Samson  (1973).  A 
series  or  papers  involving  both  practical  and  theoretical  studies  has  been 
wrift^i  by  Samson  (1973,  1977)  and  Samson  and  Olson  (1980,  1981).  The 
coherency  matrix  defined  in  the  time  domain  was  first  proposed  by  Vidale 
(1986). 

All  the  techniques  use  a  moving  window  in  either  the  frequency  or  the 
time  domain  to  define  the  polarization  characteristics.  The  moving  window 
faces  two  difficulties:  1)  when  the  polarization  is  processed  in  the  frequency 
domain,  different  arrivals  separated  in  time,  but  with  the  same  frequency 
content  will  be  mixed  together  such  that  they  can  not  be  separated  in  the  fre¬ 
quency  domain;  2)  when  polarization  is  processed  in  the  time  domain,  sig¬ 
nals  with  different  frequency  content  arriving  in  the  same  time  segment  can 
not  be  separated  in  the  time  domain. 

Jurkevics  (1988)  developed  a  rigorous  time  and  frequency  decomposition 
technique  in  which  a  50%  overlap  cosine  window  and  bandpass  filters  were 
used  in  combination.  His  technique  may  be  modified  to  overcome  the  two 
problems  mentioned  above. 

Because  of  the  complexity  of  analyzing  3-D  particle  motion,  we  will  focus 
on  the  simpler  2-D  problem  for  which  the  data  consist  of  vertical  and  inline 
time  histories  at  each  receiver  location.  In  the  absence  of  lateral  variations  in 
earth  structure,  we  would  expect  the  ground  roll  to  consist  solely  of  Rayleigh 
waves  on  these  components.  This  is  not  to  say  that  the  full  3-D  problem  is 
not  important  or  will  not  provide  interesting  results,  but  rather  that  our  pur¬ 
pose  is  to  remove  the  ground  roll  from  the  component  that  usually  has  good 
reflections,  the  vertical  component. 

THEORY 


The  purpose  of  polarization  analysis  is  to  characterize  the  phase  differ¬ 
ence  between  two  signals.  If  the  signals  differ  in  phase  by  a  constant  0  or  -it 
radians,  then  their  crosscorrelation  at  zero  lag  will  be  a  maximum.  If  there  is  a 

TT  3tt 

constant  phase  difference  of  —  or  radians,  then  their  crosscorrelation  will 

be  zero.  In  the  frequency  domain,  by  virtue  of  the  crosscorrelation  theorem, 
the  Fourier  transforms  of  the  two  crosscorrelations  will  be  pure  real  or  pure 
imaginary,  respectively,  for  the  two  cases.  On  the  other  hand,  if  a  crosscorre¬ 
lation  is  formed  between  one  signal  and  the  Hilbert  transform  of  the  other, 
then  a  maximum  crosscorrelation  occurs  when  the  original  signals  differed  in 


TT 

phase  by  —  radians.  Thus  the  formation  of  crosscorrelations  can  characterize 
the  phase  differences  between  two  signals. 

We  define  the  spectral  matrix  from  the  observed  data  as 


J  = 


<zz *  > 
<rz*  > 


<zr  > 
<rr*  > 


Jzz 

Jzr 

Jrz 

Jrr 

(1) 


where  z,  r  are  the  vertical  and  inline  seismic  field  components.  The  symbol 
<>  indicates  a  frequency  or  time  average  with  the  asterisks  denoting  the 
complex  conjugate. 
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The  input  data  are  expressed  in  complex  form  such  that  the  matrix  J  is  a 
non-negative  Hermitian,  and  contains  all  the  information  needed  to  character¬ 
ize  the  polarization  parameters.  Since  the  J  matrix  can  be  expressed  in  the 
frequency  domain  (Samson,  1973)  or  in  the  time  domain  (Vidale,  1986),  it  is, 
for  simplicity,  called  the  spectral  matrix  when  it  is  defined  as  in  (1).  To  con¬ 
struct  the  J  matrix  in  the  time  domain,  the  observed  real  time  series  is 
transformed  into  a  complex  time  series  through  the  use  of  a  Hilbert 
transform. 


For  a  quasi-monochromatic  signal  produced  by  a  physical  source,  the 
ellipticity  and  orientation  of  the  polarization  are  independent  of  time.  In 
general,  such  a  wave  field  is  partially  polarized  (Goodman,  1985).  In  particu¬ 
lar,  if  the  wave  is  a  single-frequency  wave,  it  can  be  expressed  as  the  sum  of 
one  wave  field  which  is  completely  polarized  and  one  which  is  completely 
unpolarized.  It  can  be  shown  that  this  representation  is  unique:  any  matrix 
can  be  uniquely  decomposed  into  the  form 


J  =  S  +  N  (2) 

where  S,  N  are  the  spectral  matrices  of  completely  polarized  and  completely 
unpolarized  waves  respectively. 

Usually  the  observed  spectral  matrix  contains  signal  and  noise  where  the 
signal  yields  the  completely  polarized  matrix  and  the  noise  yields  the  com¬ 
pletely  unpolarized  matrix.  If  noise  is  completely  unpolarized  and  has  the 
same  variance  on  the  two  components,  then  the  N  matrix  can  be  expressed 
by  a  diagonal  matrix  with  two  equal  elements: 


Equal  off-diagonal  elements  requires  that  the  noise  be  uncorrelated 
between  components,  and  the  two  equal  diagonal  elements  requires  that  the 
noise  energy  is  the  same  at  each  component.  In  other  words,  to  satisfy  the 
completely  unpolarized  condition  the  noise  must  be  random  and  isotropic.  A 
detailed  proof  for  the  general  case  can  be  found  in  Mott  (1986),  and,  for  the 
seismic  wave  case,  in  Shieh  (1988).  The  assumption  of  random  and  isotropic 
noise  may  not  be  true  for  some  cases  or  for  some  frequencies  ranges,  but  we 
find  that  this  assumption  works  well  for  the  seismic  data  we  examined.  If 
this  assumption  is  not  correct,  the  J  matrix  cannot  be  algebraically  decom¬ 
posed  into  the  S  and  N  components.  In  such  a  case,  noise  prior  to  the  signal 
arrival  may  be  used  together  with  the  technique  of  Samson  (1983)  to  reduce 
the  effect  of  this  noise  on  the  signal. 

We  assume  that  the  /  matrix  can  be  written  as 


J  = 


n  0 
0  n 


(4) 


For  a  completely  polarized  signal,  the  matrix  S  has  a  zero  determinant.  This 
is  easily  seen  from  the  frequency  domain  representation  of  (1)  because  in  the 
absence  of  noise  zz  rr  —  zr  z  r  =  0.  Applying  this  condition  leads  to  a 
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characteristic  eigenvalue  equation 

( Jzz  -n)(Jn  ~n)-JnJzr  =0. 

Bom  and  Wolf  (1964)  and  Mott  (1986)  showed  that  the  elements  of  the  S  and 


N  can  be  calculated  from  the  following  equations: 

i_ 

Szz  =  j(Jz!-M+j((J*+lr,?-4DETl )  ])2  (5.1) 

=  -W+yWzz  +!„f-4DETl J  ])2  (5.2) 

n  =  |(/a+/„)-{(aa  +;rr)2-4DET(  J  ])2  (5.3) 

=  U,  (5.4) 

$rz  =  Jrz  (5-5) 

where 

DEr[Jl=U  -JzrJr z 


In  equations  5.1  and  5.2,  and  are,  respectively,  the  noise-free 
smoothed  spectrum  squared  at  each  frequency,  or  energy  at  each  time,  for  the 
vertical  and  inline  components.  Equation  5.3  yields  the  smoothed  noise  spec¬ 
trum  squared  at  each  frequency,  or  energy  at  each  time  sample.  The  square 
root  of  these  quantities  represents  the  smoothed  Fourier  amplitude  spectra  in 
the  frequency  domain  or  an  RMS  amplitude  in  the  time  domain. 

Once  the  signal  and  noise  characteristics  are  isolated,  the  degree  of 
polarization  P2  is  defined  by 

p2  _  traces  _  L  4DET[  I  ]  l**7 
trace  J  {  (Jzz+Jrrf  J 

We  note  that 

0  <  P2  ^  1 


Since  DET  S=0,  one  of  the  eigenvalues  if  S  must  equal  zero,  and  the 
other  equals  the  trace  S  =  sK  +srr.  Note  that  a  physical  interpretation  of  su  is 
that  it  is  proportional  to  the  energy  of  the  coherent  signal  on  the  z- 
component,  with  a  similar  interpretation  for  sn  and  n.  Thus  P2  equals  the 
ratio  of  the  signal  energy  to  the  total  energy.  The  complex  eigenvector 
corresponding  to  the  non-zero  eigenvalue  c c  S  is 


x  = 


1/  srr/szr  J 
*1,  *2^*] 


(7) 


where  x1  and  x2  are  real  and  ij/  is  the  phase. 
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In  order  to  define  the  elliptidty,  e,  it  is  convenient  to  form  another  vec¬ 
tor 

vi  +  /v2  =  xe'4*, 

by  multiplying  the  vector  in  equation  7  by  a  constant  phase  term.  Because  of 
the  way  that  the  J  matrix  and  consequently  the  S  matrix  are  defined,  this  new 
vector  will  also  be  an  eigenvector  of  the  S  matrix.  To  define  the  elliptidty, 
Samson  and  Olson  (1980)  show  that  using 

> 

x\  sin2t|* 
x\  +  x22cos2iji 

leads  to  the  real  vectors  vj  and  v2,  with  the  constraint  v/v2=0.  The  elliptidty 
e  is  defined  as  the  ratio  of  the  Euclidian  noims  I  v  I  e  of  the  v  vectors.  The 
ratio  is  defined  to  be  less  than  1.  If  for  example,  I  vj  I  e  >  I  v2 1  e,  then 

e  =  I  v2 1  e  /  I  Vj  I  e .  (8) 

We  now  have  a  number  of  parameters  that  characterize  the  polarization 
of  the  signal,  srr,  s^,  szr,  e,  P2,  vj,  and  v2.  Any  combination  of  these  can  be 
used  to  modify  the  input  vertical  and  in-line  signals.  For  example,  the  motion 
could  be  projected  onto  the  major  axis,  P-wave  linear  motion  could  be 
rejected,  or  SV-linear  motion  could  be  rejected. 

To  construct  simple  polarization  filter  functions,  we  combine  the  two 
parameters  obtained  from  the  spectral  matrix,  which  are  the  degree  of  polari¬ 
zation  and  the  elliptidty.  A  filter  function  to  reject  noise  and  to  pass  linear 
motion  is  of  the  form 

Gt  =  Pm(l-e)n.  (9) 

A  second  filter  to  reject  linear  motion  and  noise  and  to  pass  elliptical  motion 
would  be  of  the  form 

Ge  =  Pmen  (10) 

where  m  and  n  are  integers.  The  use  of  the  term  filter  is  appropriate  since  we 
will  decompose  the  two-component  signal  into  frequency  and  time  com¬ 
ponents,  apply  the  filter  function  to  each  component,  and  then  combine  the 
results  to  yield  the  output  traces. 

PROCESSING  ALGORITHM 

The  theory  described  above  is  applicable  to  the  analysis  of  a  single 
coherent  signal  in  the  presence  of  isotropic  noise.  As  mentioned  in  the  intro¬ 
duction,  real  data  will  have  multiple  arrivals,  each  with  different  polarization 
and  frequency  characteristics.  The  algorithm  presented  here  lakes  this  into 
account.  Numbers  within  parentheses  refer  to  specific  equations  in  the  text. 


a)  Form  Complex  Z  and  R  traces  using  Hilbert  transform 

b)  Loop  over  bandpass  filter  center  frequendes 
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c)  Zerophase  bandpass  filter  complex  traces 

d)  Apply  moving  cos2  window  specified  by  filter  frequency 

e)  Construct  J  matrix  (1) 

f)  Construct  signal  matrix  (5) 

g)  Construct  degree  of  polarization  (6) 

h)  Calculate  ellipticity  (7) 

i)  Calculate  rectilinearity  (1-e) 

j)  Compute  filter  constant  Gt 

k)  If  reject  low  frequency  P  and  real  JZT  >0 

l)  G/  =  0.02  Pm(l-e)n 

m)  Else  if  reject  low  frequency  SV  and  real  Jzr<  0 

n)  Gi  =  0.02  Pm(l-e)n 

o)  Else 

p)  Gj  =  Pm(l—e)n 

q)  Multiply  windowed  segment  by  G/ 

r)  Sum  windowed  segments 

s)  Advance  window 

t)  Adjust  level 

u)  Sum  filtered,  bandpassed  traces 

v)  Change  filter  frequency 

w)  Output  polarization  filtered  inline  and  vertical  traces 


To  deal  with  multiarrival  signals,  the  vertical,  Z,  and  inline,  R,  traces  are 
decomposed  into  a  series  of  bandpass  filtered  traces  (Steps  b  and  c).  Next 
polarization  analysis  is  performed  on  segments  of  the  filtered  trace  obtained 
by  using  overlapping  cos2  windows  (Step  d).  A  cos2  window  is  used  rather 
than  the  cos  window  of  Jurkevics  (1988)  since  overlapping  windows,  shifted 
by  a  half  window  length,  will  yield  a  temporally  flat  overall  window  except  at 
the  first  and  last  segments  because  of  the  trigonometric  identity 
sin2*  +  cos2*  =  1.  The  length  of  the  cos2  window  trades  off  with  the  ability 
to  remove  incoherent  noise  (equation  5.3)  since  time  averaging  over  a  long 
window  will  reduce  the  noise,  so  that  the  extra  analysis  in  equations  5. 1-5. 3 
may  not  be  significant.  Steps  e  through  q  perform  the  polarization  analysis 
and  define  the  filter  function,  Gj.  If  this  filter  function  were  forced  to  be 
G;  =  1.0,  then  we  would  like  the  original  signal  to  pass  through  unchanged. 
Unfortunately,  since  the  bandpass  filters  are  implemented  by  successive  zero 
phase  lowpass  and  highpass  filtering  using  Butterworth  recursive  digital 
filters,  the  frequency  domain  representation  will  have  other  than  unit  gain, 
and  in  addition,  as  the  various  frequency  windows  are  summed,  the  response 
function  will  be  other  than  flat.  To  account  for  this,  we  initialize  the  process 
by  determining  the  frequency  response  of  the  filtering  operation,  but  having 
the  polarization  section  pass  the  input  signal  completely,  and  compute  a  nor¬ 
malizing  gain  factor.  This  normalization  is  applied  in  Step  t.  The  result  of 
this  is  that  a  linearly  polarized  chirp  or  single  frequency  signal  is  passed 
without  any  change  in  amplitude  or  phase. 
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Steps  e  through  i  follow  the  equations  in  the  text.  However,  one  other 
feature  is  built  into  the  algorithm,  that  makes  use  of  the  fact  that  vector  parti¬ 
cle  motion  data  are  acquired.  If  vertical  ground  motion  upward  and  inline 
ground  motion  away  from  the  source  are  defined  positive,  then  the  sign  of 
the  real  part  of  Jzr  can  be  used  to  specify  whether  the  ground  motion  is  of  P 
or  SV  character,  according  to  whether  the  sign  is  positive  or  negative,  respec¬ 
tively.  In  Steps  k  through  p  we  can  specify  a  rejection  of  P  or  SV  motion  at 
the  same  time  as  we  reject  incoherent  noise  and  non-linear  particle  motion. 
This  is  applied  to  remove  low  frequency  linearly  polarized  ground  roll  com¬ 
ponents  in  the  signal. 

The  parameters  used  for  default  processing  are  m  =  2,  n  =  6,  signal  pro¬ 
cessing  with  center  frequencies,  fc ,  up  to  100  Hz,  and  a  bandwidth  of  10A / , 
where  A /  is  the  Fast  Fourier  transform  sampling  frequency,  obtained  as  a 
result  of  performing  the  Hilbert  transform  through  frequency  domain  opera¬ 
tions.  The  choice  of  m  and  tt  is  not  critical;  the  values  used  were  chosen 
after  some  experimentation  with  real  data  sets.  The  time  domain  moving  cos2 
window  has  a  window  width  equal  to  2  /  fc .  If  rejection  of  low  frequency  P 
or  SV  components  is  requested,  it  is  done  for  frequencies  less  than 
fL  =  35Hz. 

DATA  PROCESSING 

Several  different  data  sets  will  be  presented  to  test  the  operation  of  the 
technique  as  well  as  to  consider  its  applicability  under  different  conditions  of 
reflection  signal  to  ground  roll  noise  level. 

Synthetic  Data  Set 

Figure  1  presents  some  synthetic  data  sets  together  with  the  results  of 
our  polarization  analysis.  The  synthetic  vertical  component  time  history  con¬ 
sisted  of  four  30  Hz  Ricker  wavelets  centered  at  two  way  travel  times  of  0.20, 
0.35,  0.50,  0.65  and  0.80  seconds.  The  synthetic  inline  time  history  consisted 
of  30  Hz  Ricker  wavelets  phase  shifted  by  0°,  45°,  90°,  135°  and  180°,  at  the 
corresponding  times.  The  vertical  component  wavelets  had  a  peak  amplitude 
of  1.0.  In  addition  zero  mean  random  noise  with  a  Gaussian  distribution  and 
a  standard  deviation  of  A,  where  A  is  an  amplitude  level,  was  added  to  the 
polarized  signals.  Since  polarization  analysis  would  usually  be  applied  to  data 
prior  to  the  NMO  stack,  it  is  also  useful  to  see  how  it  affects  the  stack. 

Figure  1  consists  of  four  groups  of  six  traces.  In  each  group,  the  sixth 
and  rightmost  trace  is  a  simple  stacked  sum  of  the  preceding  five  traces.  Fig¬ 
ure  la  corresponds  to  the  five  vertical  component  time  history  realizations  for 
A  =0.2  and  the  resulting  stack.  Figure  lb  shows  the  result  of  polarization 
filtering  each  of  the  five  input  traces  and  the  corresponding  stacked  output. 
We  see  that  the  polarization  filter  did  reject  the  wavelets  with  phases  other 
than  0°  or  180°.  In  addition,  there  is  some  improvement  in  the  signal  to  noise 
ratio  within  the  individual  traces,  indicating  some  value  in  attempting  to 
define  the  isotropic  noise  component  in  the  signal  .  The  stack  of  the  polariza¬ 
tion  filtered  traces  has  slightly  less  noise  than  the  stack  of  the  input  traces. 
Figures  1c  and  Id  compare  the  input  and  polarization  filter  output  traces  and 
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stacks  for  the  case  when  the  noise  is  doubled,  A  =0.4.  The  stacked  traces 
again  show  the  benefit  of  noise  reduction  by  stacking,  but  again  the  non- 
linearly  polarized  wavelets  are  rejected  by  the  filter.  Another  simulation  was 
run  with  A  =0.8,  but  the  results  indicated  a  reduction  in  the  ability  to  reduce 
the  noise  and  also  to  reduce  the  non-linearly  polarized  signal. 

Oklahoma  Field  Data 

As  part  of  an  sensor  evaluation  test,  four  data  sets  were  collected  at  a 
field  test  site  in  Oklahoma.  Both  dynamite  and  vertical  vibrator  sources  were 
used  with  receiver  group  arrays  spaced  15.2  meters  apart  between  152.4  and 
441.9  meters  from  the  source.  Each  group  array  consisted  of  six  geophones. 
The  vibrator  acted  at  the  surface  while  the  dynamite  was  buried  at  a  depth  of 
15.2  meters.  The  four  data  sets  acquired  are  as  follows: 

DynGrp  -  dynamite  source  with  a  30.5  meter  geophone  group  array 
DynNoGrp  -  dynamite  source  with  a  1  meter  geophone  group  array 
VibGrp  -  vibrator  source  with  a  30.5  meter  geophone  group  array 
VibNoGrp  -  vibrator  source  with  a  1  meter  geophone  group  array 

The  receiver  group  arrays  were  not  designed  to  maximally  reduce  the  surface 
waves;  at  most  they  reduce  the  surface-wave  signal  by  6-10  db  over  the 
usable  signal  range  (5-85  Hz). 

Figure  2  presents  the  acquired  vertical  component  time  histories,  which 
consist  of  1000  ms  of  data  sampled  at  4  ms  intervals.  Within  each  panel,  the 
trace  on  the  left  is  at  a  distance  of  152.4  m  and  the  one  at  the  right  is  at  a  dis¬ 
tance  of  441.9  m.  An  automatic  gain  correction  with  a  500  ms  window  has 
been  applied  to  each  trace.  A  well  developed  set  of  reflections  is  seen  at  a 
two  way  traveltime  of  0.550  seconds.  The  sections  are  ordered  so  that  the 
quality  of  the  reflection  decreases  from  left  to  right. 

Figure  3  presents  a  comparison  of  the  input  and  output  traces  for  the 
DynGrp  data  set.  Figures  3a  and  3b  are  the  input  vertical  and  inline  data 
sets,  respectively,  while  Figures  3c  and  3d  are  the  filtered  vertical  and  inline 
data  sets,  respectively.  The  reflection  at  0.55  seconds  has  been  enhanced  on 
the  vertical  components.  In  addition,  much  of  the  reverberation  following  the 
direct  arrival  has  been  associated  with  non-linear  particle  motion  and  subse¬ 
quently  removed.  Some  low  frequency  ground  roll  is  still  present  in  Figure 
3c.  The  inline  output.  Figure  3d,  does  not  show  any  good  reflections,  which 
it  should  not  for  rays  nearly  vertically  incident  at  the  free  surface.  The  low 
frequency  ground  roll  and  the  reverberations  following  the  direct  P  arrival 
have  been  reduced,  though.  Because  of  the  lack  of  any  distinct  arrivals  on 
the  inline  component,  we  will  not  present  further  displays  of  this  component 
with  this  data  set. 

Figure  4  presents  the  results  for  the  vertical  component  of  the  data  set 
DynNoGrp.  In  this  case  Figure  4a  represents  the  input  vertical  component 
data  set.  Figure  4b  represents  the  polarization  filter  output.  Figure  4c  is  the 
polarization  filter  output  with  the  additional  rejection  of  low  frequency  P- 
wave  motion,  and  Figure  4d  is  the  polarization  output  with  the  rejection  of 
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the  low  frequency  SV-wave  motion.  The  rejection  of  low  frequency  P-wave 
motion  on  the  vertical  component.  Figure  3c,  is  quite  effective.  The  data  pro¬ 
cessing  resolves  the  nature  of  the  reflection  at  0.55  seconds  very  well. 

Finally,  Figure  5  presents  the  polarization  filtered  vertical  component 
traces  for  each  of  the  above  data  sets,  arranged  in  order  of  increasing  surface 
wave  signal  in  the  original  data  set.  In  addition  all  low  frequency  P-wave 
motion  at  frequencies  less  than  35  Hz  was  attenuated.  Figures  5a,  5b,  5c  and 
5d  correspond  to  the  DynGrp,  DynNoGrp,  VibGrp,  and  VibNoGrp  data  sets 
respectively.  In  all  cases,  there  is  a  definite  reduction  in  the  ground  roll, 
together  with  corresponding  relative  enhancement  of  the  reflections.  From 
this  we  conclude  that  the  technique  works  as  designed,  and  that  the  algo¬ 
rithm  does  not  introduce  any  spurious  arrivals. 

Figure  6  presents  another  Oklahoma  data  set,  DYN.  These  data  were  col¬ 
lected  from  a  4  kg  dynamite  source  buried  at  a  depth  of  16  m  and  were 
recorded  at  distances  between  20  and  1000  meters.  Each  trace  is  due  to  a  sin¬ 
gle  geophone,  so  that  the  group  array  is  of  length  0  m.  Figure  6a  gives  the 
input  vertical  component  time  histories  and  Figure  6b  gives  the  polarization 
filter  output,  with  the  low  frequency  P-wave  rejected.  It  is  interesting  to  note 
that  many  of  the  reverberations  apparent  in  the  initial  record  are  significantly 
reduced  in  amplitude.  This  would  lead  to  improved  results  during  NMO 
velocity  analysis. 

Permafrost  Field  Data 

Data  sets  in  the  arctic  regions  are  notorious  for  poor  signal  to  noise,  and 
often  are  excellent  candidates  for  phase  match  filters  (Barton  et  al.,  1986; 
McConnell  et  al.,  1986;  and  Beresford-Smith  and  Rango,  1988)  Our  data  set, 
the  vertical  component  traces  are  shown  in  Figure  7a,  was  acquired  on  land, 
and  is  not  overwhelmed  by  the  dispersive  flexural  waves  typical  of  data 
acquisition  on  winter  ice  sheets.  The  traces  range  in  distance  from  1542  m  on 
the  left  to  100  m  on  the  right.  There  are  some  very  strong  undispersed 
arrivals  with  group  velocities  of  1400  m/s  and  300  m/s.  The  large  amplitude 
of  these  arrivals  is  indicated  by  the  apparent  muting  introduced  immediately 
above  them  due  to  the  AGC  process.  These  signals  are  also  very  narrow 
band  with  a  center  frequency  near  15  Hz.  The  results  of  processing  the  verti¬ 
cal  component  and  not  rejecting  either  the  low  frequency  P  or  SV  arrivals  is 
shown  in  Figure  7b. 

Comparing  the  polarization  filter  output  to  the  input  time  history,  we 
note  that  the  direct  arrival  is  eliminated,  that  the  overall  record  appears 
simpler  due  to  the  elimination  of  a  what  may  be  reverberations.  The  prom¬ 
inent  air  wave  is  not  reduced  in  amplitude,  indicating  that  in  this  case  its  par¬ 
ticle  motion  is  highly  linear. 

CONCLUSIONS 

A  polarization  filter  technique  has  been  presented  for  use  with  vertical 
and  inline  data  traces  for  the  purpose  of  reducing  ground  roll.  The  technique 
works  well  when  the  ground  roll  components  are  associated  with  non-linear 
particle  motion,  but  does  pass  linearly  polarized  components  of  the  ground 
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roll.  This  is  not  unexpected  since  the  ground  roll  elliptidty  can  vary  signifi¬ 
cantly  from  near  circular  polarization  to  almost  linear  polarization  as  a  func¬ 
tion  of  frequency  and  mode. 

The  polarization  technique  also  reduced  the  reverberation  following  the 
direct  P-wave  arrival  at  large  offset,  and  thus  may  be  useful  in  enhancing 
other  non-vertically  incident  arrivals  that  are  important  for  imaging  subsur¬ 
face  structure. 

The  theory  of  polarization  filtering  was  developed  only  sufficiently  to 
accomplish  the  task  of  2-D  filtering.  Other  information,  such  at  the  direction 
of  linear  particle  motion  at  the  surface,  a  function  of  the  angle  of  incidence  at 
the  free  surface,  was  not  used.  This  other  information  may  be  useful  for 
advanced  imaging  of  the  vector  wavefield,  especially  in  the  case  of  three  com¬ 
ponent  data  sets,  for  which  it  may  be  possible  to  focus  the  data  set  upon  a 
specific  set  of  offline  reflectors. 

The  improvement  in  signal  to  noise  ratio  is  not  as  profound  with  the  sur¬ 
face  sources.  On  the  other  hand,  in  the  case  of  initially  good  reflection  signal 
to  ground  roll  noise,  as  is  the  case  for  buried  dynamite  sources,  the  polariza¬ 
tion  filter  technique  yields  very  good  results  even  with  minimal  group  arrays 
(Figure  6).  Thus  the  full  use  of  the  vectorial  nature  of  ground  motion  may 
lead  to  less  intensive  data  acquisition  strategies  with  no  noticeable  decrease  is 
processed  line  quality. 
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Fig.  2.  Input  data  sets  from  the  Oklahoma  test  site,  a)  dynamite  source,  group 
array  DynGrp;  b)  dynamite  source  1  m  group  array  DynNoGrp;  c) 
vibrator  source,  group  array  VibGrp;  d)  vibrator  source,  1  m  group 
array  VibNoGrp.  All  traces  are  AGC'd  with  a  500  ms  window. 
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.  Polarization  processing  of  the  DynGrp  data  set.  All  traces  in  the 
display  have  had  an  AGC  applied  with  a  500  ms  window,  a)  Input 
vertical  component  time  histories,  b)  Input  inline  component  time  his¬ 
tories.  c)  Output  vertical  component  time  history,  d)  Output  inline 
component  time  history. 
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Fig.  6.  Polarization  processing  of  the  DYN  data  set.  All  traces  in  the  display 
have  had  an  AGC  applied  with  a  500  ms  window.  In  addition,  all 
have  been  processed  to  reject  P-wave  motion  at  frequencies  less  than 
35  Hz.  a)  Input  vertical  component  time  history,  b)  Output  vertical 
component  time  history. 
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.  Vertical  component  time  histories  for  the  permafrost  data  set.  All 
traces  in  the  display  have  had  an  AGC  applied  with  a  500  ms  window, 
a)  Input  time  histories,  b)  Output  time  histories. 
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